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a 


* 


a 


a 


asc*otbcc»af  cc 


8 

8 

8* 

Y 

Yi 

Yt 

YiJ 

6ij 

6ij 

<5 

eU 


constant;  value  of  the  Eshelby  S-tensor 

relative  vertical  displacement  between  the  centers  of  two 
spheres  in  contact 

angle  of  orientation  of  major  principal  stress 

length  of  edge  of  representative  volume  of  the  sc,  bcc 
and  fee  arrays  respectively 

constant  (=  dT/dN) 

angle  of  orientation  of  the  element 
constant;  value  of  the  Eshelby  S-tensor 
engineering  shear  strain 

average  shear  strain  experienced  by  inclusion  i 
threshold  strain 

engineering  shear  strains  (components  of  the  strain 
t  ensor) 

Kronecker  delta 
displacement 

horizontal  displacement  between  the  centers  of  two 
spheres  in  contact 

strain  tensor 


e  p 

ei.ejj  :  elastic  and  plastic  strains  respectively 

ev  ■  average  volumetric  strain  experienced  by  inclusion  i 


ev 


0 

®z 


:  average  macroscopic  volumetric  strain 
:  constant  (=f/8) 

:  angle  of  Incidence  of  S-wave 


i*a  I4*  4*«  * 1.1  l,» J.l  i.) 


_  ,  ,  ......  ,  ...  ...  ___  ...  ...  l|«  ft*  ■ 


coefficient  of  proportionality  in  plastic  flow  rule 


constant,  parameter  of  the  probability  density  function 
of  the  void  ratio,  p(e) 


Poisson's  ratio 


macroscopic  Poisson's  ratio 


Poisson's  ratio  of  the  spheres 


3.14159 


mass  density 


Oil  (1=1, 2, 3)  :  normal  stress  applied  to  regular  array 


stress  tensor 


tflj  (i  ,j  =  l ,2,3) :  shear  stress  applied  to  regular  array 


applied  (macroscopic)  state  of  stress 


axial  stress 


°a»  °b 


principal  stresses  in  directions  of  propagation  and 
polarization 


principal  stress  in  direction  of  no  propagation/no 
polarization 


isotropic  confining  pressure 


mean  effective  stress 


oi  (1=1, 2, 3)  :  principal  stress 


(1=1, 2, 3)  :  effective  principal  stress 


isotropic  stress  experienced  by  inclusion  i 


applied  isotropic  stress 


applied  shear  stress 


shear  stress 


t*  :  shear  stress  at  failure  (=f<j0) 

av,on  :  vertical  and  horizontal  stresses 


Latin  Letters 


A 

a 

a 

C 

cijM, 

^n»Ct 

cn 

ci 

D,  Di:j 
D,D 

D10 
d^D  ,dD 
E 

Es 

e 

erain»emax 

e 

F 

£(^ij) 

t 

FEM 


:  area  of  face  of  elementary  volume 
:  constant  (0.15  _<  a  <_  0.23) 

:  radius  of  area  of  contact  between  two  spheres  in  contact 
:  constant 
:  Compliance  Matrix 

:  Normal  and  Tangential  Compliance  respectively 
:  Uniformity  Coefficient 
:  volume  concentration 
:  Constrained  Modulus 

:  Displacement  between  centers  of  two  spheres  in  contact 
:  percent  passing 
:  Displacement  increment 
:  Young's  Modulus 
:  Young's  Modulus  of  the  spheres 
:  void  ratio 

:  minimum  and  maximum  void  ratio  respectively 
:  mean  void  ratio 
:  constant 
:  yield  function 
:  friction  coefficient 


:  finite  element  method 


:  shear  modulus  of  Inclusion  i  (i  =  1,...,N) 


„-v  ’.u -s-  \-v  V  »■».’  V-  17 


Gi 


G* 


G 

GS 

G 

Gmax 

g(aij) 

W.Ho.Hp 

M 

*i 

K* 

K 

k 

L 

L* 

M 

ra 

ma.mb.mc 


n 

nmln»  nmax 
N 

N 

N 


macroscopic  shear  modulus 

secant  shear  modulus 

shear  modulus  of  the  spheres 

normalized  shear  modulus 

shear  modulus  at  very  small  strains 

plastic  potential 

tangential  elastoplastic  and  elastic  moduli 

first  invariant  of  the  stress  tensor 

Bulk  Modulus  of  inclusion  i  (i  =  1,...,N) 

macroscopic  Bulk  modulus 

stress  ratio  (K  =  oz/a0) 

constant  (k  *  (2-vs)/2(l  vs)) 

constant  (=  T/fN0) 

constant  (=  T*/fNQ) 

constant 

constant 

constants 

constant 

porosity 

mean  porosity,  macroscopic  porosity 

minimum  and  maximum  porosity  respectively 

Number  of  cycles 

Number  of  phases  or  materials 

Normal  force 

xx  i 


4 


V 

I 

3 


3 

v 

a 

Ss 


WjJ 

Wjl 


Vj 

V 

S 


vf 


1 


Contact  forces  at  local  coordinate  system 


Contact  forces 


Normal  contact  force  due  to  application  of  o0 


constant 


contact  normal  force 


constant 


contact  force 


atmospheric  pressure 
applied  forces 

probability  density  function  of  x 
radius  of  spheres 


elastic  constants 


Stiffness  matrix 


Tangential  contact  forces  in  x  and  y  directions 


dTn,dTt 


outward  normal  and  tangential  components  to  the  yield 
surface  of  the  applied  tangential  force  increment 
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coordinate 


tangential  force 


tangential  force  at  failure  (=  fN0) 


P  and  S  wave  velocities 


rod  wave  velocity 
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ABSTRACT 


The  need  for  a  micromechanical  approach  to  modeling  the  stress- 
strain  response  of  granular  soil  is  discussed  and  justified.  This  work 
focuses  on  the  small  shear  strain  (y  _<  0.01%)  behavior,  and  investigates 
the  validity  of  modeling  analytically  uniform,  rounded-grained  quartz 
sand  by  arrays  of  identical  elastic,  rough,  quartz  spheres.  As  a  first 
step,  the  stress-strain  properties  of  six  regular  arrays  of  spheres  are 
studied  in  some  detail,  with  focus  on  isotropic  and  and  biaxial  boundary 
loading. 

An  analytical  procedure  is  established  for  determining  the  elastic 
moduli  of  a  random  assemblage  of  equal,  elastic,  rough  spheres  of 
arbitrary  mean  porosity,  subjected  to  isotropic  confining  pressure.  The 
procedure  uses  the  properties  of  the  regular  arrays  already  described,  it 
accounts  for  the  spatial  distribution  of  porosity,  and  it  calculates  the 
macroscopic  moduli  through  the  Self  Consistent  Method.  The  procedure  was 
applied  to  compute  the  shear  and  bulk  moduli  of  assemblages  of  quartz 
spheres,  which  were  then  compared  to  static  and  dynamic  measurements  on 
uncycled  and  heavily  precycled  quartz  sands  reported  in  the  literature. 
Although  the  theoretical  sands  are  significantly  stiffer  than  the  actual 
soils,  excellent  agreement  was  found  with  resonant  column  measurements  on 
heavily  precycled  Ottawa  sand  at  small  strains. 

Finally,  a  new  two-dimensional  model  of  the  stress-strain  behavior 
of  granular  soil  at  small  strains  is  presented.  The  model  is  based  on  an 
incremental  solution  to  the  contact  problem  of  two  equal,  elastic,  rough 
spheres  and  is  implemented  through  nonlinear  finite  element  techniques. 
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The  results  of  numerical  experiments  conducted  on  this  idealized 
aggregate  are  compared  to  laboratory  data  on  the  static  and  cyclic  small 
strain  behavior  of  actual  sand,  as  well  as  to  recent  compressional  wave 
velocity  measurements  on  anisotropically  consolidated  dry  sands,  with 
good  agreement.  These  measurements,  performed  at  the  large  cubic 
triaxial  facility  at  the  University  of  Texas,  have  shown  that  the  P-wave 
velocity  depends  only  on  the  principal  stress  parallel  to  the  direction 
of  wave  propagation;  this  finding  was  also  predicted  by  the  simulation. 

Concluding,  the  hypothesis  that  certain  aspects  of  the  behavior  of 
granular  soil  are  due  to  the  particulate  nature  of  the  soil  is  justified. 
These  aspects  cannot  be  interpreted  and  reproduced  analytically  unless 
this  particulate  nature  is  taken  into  account. 


CHAPTER  1 


INTRODUCTION 

The  main  objective  of  this  work  is  to  present  a  simple,  yet 
rigorous,  particulate  mechanics  model  of  the  stress-strain  response  of 
granular  soil  under  small  shear  strains,  y»  of  the  order  of  10 “2%  or 
less.  The  proposed  model  idealizes  sand  as  a  combination  of  regular 
arrays  of  elastic,  rough  spheres  and  uses  Mindlin’s  formulation  for  the 
contacts . 

"Elastic  constants"  of  interest  at  very  small  strains^*)  include  the 
shear  and  bulk  moduli  and  the  Poisson’s  Ratio(s).  Experimental  results 
and  basic  considerations  indicate  that  these  "constants"  depend  on  both 
the  void  ratio  of  the  soil  and  the  state  of  confining  stresses.  The 
variations  of  these  moduli  and  of  the  damping  of  the  soil  with  an  applied 
shear  strain  up  to  the  threshold  are  also  of  interest,  as  is  the  value  of 
the  threshold  strain  itself  at  which  gross  sliding  occurs  at  the  grains’ 
contacts . 

These  small  strain  soil  parameters  are  very  important  in  geotechni¬ 
cal  engineering  problems  involving  cyclic  loading  or  wave  propagation  in 
the  soil,  such  as:  ocean  wave  loading,  soil  structure  interaction,  site 
response,  ground  settlement  and  liquefaction  during  earthquakes.  Due  to 
this,  a  great  number  of  experimental  studies  of  small  strain  behavior 


The  concepts  of  small  and  large  shear  strains  as  used  here  are 
consistent  with  usual  Soil  Dynamics  terminology,  but  they  do  not 
coincide  with  that  of  traditional  Soil  Mechanics,  where  much  greater 
strains  of  about  1%  or  larger  are  usually  of  interest. 
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have  been  performed,  and  correlations  have  been  developed  for  practical 
use.  Especially  important  are  the  equations  for  the  shear  modulus  at  very 
small  strains,  Gmax,  in  sands  developed  by  Hardin  and  Richart  (1963)  and 
Seed  and  Idriss  (1970),  be  used  on  the  assumption  that  these  soils  can  be 
treated  as  elastic  isotropic  solids.  In  both  correlations,  discussed  in 
more  detail  in  Section  2,  G^x  =  with  aj,  02»  a 3  being  the 

effective  principal  stresses,  aQ  =  (ai  +  02  +  03)/3  is  the  mean  effective 
stress  and  A  is  a  soil  constant  which  depends  on  void  ratio  or  relative 
density.  Both  correlations  assume  that  G^x  (and  thus,  also,  the  shear 
wave  velocity,  Vs  »  (Gmax/p)^^)»  is  the  same  for  isotropically  or  aniso- 
tropically  loaded  sand,  provided  that  the  mean  stress  aQ  is  the  same; 
also,  both  correlations  assume  that  for  the  anisotropic  loading  case  Gmax 
and  Vs  do  not  change  with  direction. 

These  assumptions  for  Gmax  in  sands  have  been  challenged  more  re¬ 
cently  by  the  experimental  results  obtained  by  Roesler  (1979),  Knox  et 
al.  (1982)  and  Yu  and  Richart  (1984),  as  discussed  in  Section  2.  There¬ 
fore,  a  main  motivation  for  this  work  was  the  need,  suggested  by  those 
experimental  findings,  for  a  fresh  approach  to  our  basic  understanding  of 
Gmax  and  other  small-strain  soil  parameters.  Some  preliminary  analytical 
results  previously  obtained  by  Dobry  et  ai.  (1982),  had  shown  that  a 
particulate  mechanics  approach  was  very  well  suited  to  this  purpose,  and 
should  be  the  basis  of  this  fresh  approach. 

The  large  strain  (0.012  >  y)  behavior  of  granular  soils  is  also  very 
important  in  engineering  problems  involving  cyclic  loading,  and  espe¬ 
cially  those  related  to  earthquakes.  At  these  strains,  the  stress-strain 


behavior  becomes  strongly  nonlinear  and  hysteretic,  and  rearrangement  of 
particles  take  place  producing  phenomena  such  as  densif ication  and  pore 
water  pressure  build-up  (Silver  and  Seed,  1971;  Youd,  1972;  Dobry  et  al., 
1982;  National  Research  Council,  1985).  A  number  of  continuum  mechanics 
models,  based  mostly  on  the  Incremental  Theory  of  Plasticity,  try  to 
simulate  this  behavior  and  have  been  proposed  for  soils,  as  discussed  in 
Chapter  3. 

A  summary  literature  review  of  previous  relevant  particulate 
mechanics  studies  is  presented  in  Chapter  4.  Many  of  these  past  inves¬ 
tigations  have  focused  on  the  very  large  strain  (y  >  1%)  and  strength 
behavior  of  granular  soils;  at  those  very  large  strains,  gross  sliding 
and  rolling  of  the  grains  are  main  contributors  to  the  overall  strain, 
while  the  elasticity  of  the  particles  and  contacts  play  a  minor  or 
negligible  role.  On  the  other  hand,  for  the  small  to  large  strain  ranges 
of  Interest  of  the  proposed  research,  the  elasticity  of  the  particles  and 
the  details  of  the  force-displacement  response  at  the  contacts  are  very 
significant  factors.  The  discussion  in  Chapter  4  includes  a  general 
model  recently  proposed  for  the  force-displacement  response  at  the 
contact  between  two  identical  elastic  spheres  (Seridi  and  Dobry,  1984, 
Dobry  et  al. ,  1987) . 

The  results  of  the  present  research  are  discussed  in  Chapters  5 
through  9.  Chapter  5  presents  a  deta. led  study  of  the  differential 
stress-strain  relations  for  various  regalar  arrays  of  spheres.  Chapter  6 
describes  an  application  of  these  findings,  using  the  Self-Consistent 
Method,  to  a  random  arrangement  of  regular  arrays  subjected  to  isotropic 
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boundary  loading,  and  with  the  arrangement  having  an  arbitrary  macro¬ 
scopic  void  ratio.  Chapter  7  presents  the  formulation  of  a  two- 
dimensional  simulation  using  nonlinear  finite  element  techniques  and  the 
contact  model  of  Chapter  4.  The  results  of  this  simulation  appear  in 
Chapters  8  and  9,  where  an  excellent  agreement  is  noted  between  the 
simulation  and  laboratory  measurements. 


CHAPTER  2 


LABORATORY  MEASUREMENTS  ON  SANDS  AT  SMALL  STRAINS 

Starting  around  1960,  a  number  of  cyclic  and  dynamic  laboratory 
measurements  have  been  performed  to  determine  the  stress-strain  behavior 
of  granular  arrays  and  of  natural  sands  at  small  strains.  Properties 
studied  have  included:  i)  maximum  shear  modulus  at  very  small  strains, 
Gmax;  ii)  the  variation  of  secant  modulus,  G,  with  shear  strain,  y;  iii) 
the  Poisson's  ratio  of  the  soil;  iv)  the  variation  of  shear  damping  ratio 
with  strain;  and  v)  the  threshold  shear  strain,  yt ,  at  which  densifica- 
tion  and  pore  pressure  buildup  start.  Many  of  these  tests  have  been  con¬ 
ducted  in  a  triaxial  cell,  on  sand  specimens  isotropically  or  anisotropi- 
cally  consolidated  under  a  biaxial  stress  state  (02  =  03  or  02  =  a^), 
with  the  small  strain  measurements  performed  using  the  pulse  method,  the 
resonant  column  or  cyclic  torsional  techniques,  and  with  particular 
emphasis  on  shear  modulus  determinations.  Important  results  and  state- 
of-the-art  summaries  of  these  modulus  measurements  have  been  presented  by 
Duffy  and  Mindlin  (1957),  Hardin  and  Richart  (1963),  Lawrence  (1965), 
Richart  et  al.  (1970),  Seed  and  Idriss  (1970),  Hardin  and  Drnevich 
(1972),  Woods  (1978),  Iwasaki  et  al.  (1978),  and  Tatsuoka  et  al.  (1979). 

As  mentioned  before,  "small  strains"  are  defined  here  by  the  condi¬ 
tion  y  <  Yt>  as  in  this  range  the  original  geometry  of  the  granular  array 
or  sand  remains  essentially  unchanged,  with  very  few  or  no  particles 
experiencing  gross  sliding  or  rolling,  and,  thus,  with  the  macroscopic 
strain  of  the  array  being  controlled  by  the  elastic  deformations  of  the 
particles  and  by  localized  slips  within  the  areas  of  contact  areas 


between  particles.  In  many  sands,  yt  a  10“^*  =  10~^%;  measurements  and 
studies  of  Yt  have  been  presented  by  Drnevich  and  Richart  (1970),  Youd 
(1972),  Pyke  (1973),  Dobry  et  al.  (1980,  1981,  1981a,  1982),  Dyvik  et  al. 
(1982),  Oner  (1984)  and  National  Research  Council  (1985). 

On  the  basis  of  laboratory  measurements,  Hardin  and  Richart  (1963) 
proposed  the  following  expression  for  Gma3C : 


f(e)(a0)0-5 


where  e  *  void  ratio,  f(e)  *  2630( 2. 1 7-e )2/  ( 1-t-e)  for  round-grained  sands, 
and  f(e)  =  1230(2. 97-e)^/ (1+e)  for  angular-grained  sands;  both  and 

aQ  are  in  psi  in  Eq.  1. 

Seed  and  Idriss  (1970)  proposed  the  alternative  expression: 


Gmax  =  1.000  K2max  (oo)0-5 


where  G^ax  and  oQ  are  in  psf,  and  K2max  Is  a  function  of  relative 
density. 

Eqs.  1-2  reflect  the  conclusion  of  these  and  other  studies,  that 
Gmax  and  Vs  are  mainly  a  function  of  void  ratio  or  relative  density,  and 
of  the  mean  effective  normal  stress  o0.  Other  variables,  such  as  static 
shear  stresses,  stress  history,  stress  path  (compression  versus  extension 
loading),  frequency  of  cyclic  loading,  degree  of  saturation,  were  found 
to  have,  either  a  small  effect  or  no  effect  at  all  (Richart  et  al.,  1970, 
Yu  and  Richart,  1984).  One  exception  is  that  a  large  number  of  shear 
ptostraioing  cycles  at  strains  larger  than  the  threshold  was  found  to 
increase  G,,,-^  significantly  (Drnevich  and  Richart,  1970). 
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It  is  useful  to  make  explicit  some  of  the  implications  of  Eqs.  1-2 
for  anisotropically  loaded  dry  sand,  either  for  the  general  "triaxial" 
case  in  which  ai  *  a 2  *  03  or,  as  very  usually  the  case  in  the  field, 
the  "biaxial"  case,  o\  2  ~  03*  (The  bars  have  now  been  dropped  from 
the  stresses,  as  for  a  dry  sand,  a  =  a).  These  implications  are: 

i)  Gmax  and  Vs  depend  equally  on  03,  02  anc*  a3* 

ii)  For  a  shear  wave  propagating  in  the  sand,  the  value  of  Vs  is  the 
same  whatever  the  directions  of  propagation  and  polarization  of 
the  wave. 

Implication  ii)  is  equivalent  to  assume  that,  at  very  small  strains, 
the  anisotropically  loaded  sand  can  be  modeled  as  an  isotropically 
elastic  material,  defined  by  the  two  elastic  constants  Gmax  and  Poisson's 
Ratio,  v.  Under  the  usual  additional  assumption  that  for  a  dry  sand,  v  = 
0.3  to  0.4  is  a  constant  independent  of  confining  stresses,  a  set  of 
implications  similar  to  i)  and  ii)  can  be  obtained,  but  now  for  a 
dilational,  P-wave,  travelling  in  this  dry  sand.  If  we  call  D,  Vp  the 
constrained  modulus  and  P-wave  velocity,  respectively,  the  corresponding 
implications  are: 

iii)  D  and  Vp  depend  equally  on  aj,  a 2  and  03,  with  D  a(o0)^’^  and  Vp 
a(o0)^*25* 

iv)  The  value  of  Vp  is  the  same  whatever  the  direction  of  propagation 
of  the  P-wave. 

Five  recent  experimental  studies  have  attempted  to  verify  in  detail 
this  formulation  by  Hardin/Richart  and  Seed/Idriss,  and  specifically  the 
validity  of  implications  i)  through  iv)  above  for  anisotropically  loaded 
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dry  sand. 

Schmertmann  (1978)  measured  Vp  and  Vg  in  several  directions  in  a 
large  dry  sand  specimen  (4  ft  diameter  by  4  ft  high).  In  these  tests,  a 
biaxial  state  of  stress  could  be  achieved,  03  =  ov  *  03  ~  03  =  with 
av,  ah  =  vertical,  horizontal  stresses,  with  the  stresses  varying  between 
5  to  20  psi,  and  with  a  stress  ratio,  01/03  =  1  to  3.  He  found  that 
there  was  a  slight  amount  of  inherent  anisotropy  (different  wave 
velocities  in  the  horizontal  and  vertical  directions  when  av  =  o^).  He 
also  found  that  for  constant  o0  =  l/3(ov  +  2  a^)  and  variable  03/03,  Vs 
varied  less  than  10X,  thus  verifying  the  basic  Hardin/Richart  assumption 
as  a  first  approximation  for  Vs  in  this  biaxial  case.  However,  Vp 
was  strongly  affected  by  03/03.  The  results  suggested  that,  for  P-waves 
propagating  in  the  vertical  direction,  Vp  depended  more  on  ov  than  on 


Roesler  (1979)  measured  Vs  using  a  1  ft^  cubical  dry  sand  sample. 
In  these  tests,  a  true  triaxial  state  of  stresses  was  achieved.  Test 
pressures  ranged  from  5.8  psi  to  23  psi,  with  03/03  =  1  to  1.8.  He 
propagated  the  shear  waves  along  either  of  the  principal  stress 
directions  (oa),  with  particle  motions  polarized  in  another  principal 
direction  (03,).  The  third  principal  direction,  or  out-of-plane 
direction,  is  neither  a  direction  of  propagation  nor  polarization  (oc). 
Roesler  found  that  his  results  followed  the  law: 


Vs  =*  B  0  0.149  0,0 .107  0 
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where  B  =•  constant.  These  results  are  Illustrated  in  Fig.  1.  For 
isotropic  confinement  (0o=<ja=ob=ac)  they  do  confirm  the  Hardin-Richart 
law  that  Gmax  a(o0)®*^  and  Vs  a(o0)®'25,  as  0.149  +  0.107  =  0.256. 
However,  for  the  general  case,  Eq.  3  contradicts  Eqs.  1-2,  in  that  now  Vs 
is  completely  independent  of  oc.  Also,  Roesler's  results  for  this  case 
indicate  that  Vs  is  a  function  of  direction  and  the  sand  cannot  be 
treated  as  an  isotropic  elastic  material;  therefore,  more  than  two 
elastic  constants  are  necessary  to  define  it. 

Stokoe  et  al.  (1980)  developed  at  the  University  of  Texas  at  Austin 
(U.T.)  a  large  scale,  7  ft^  cubical  triaxial  facility,  for  the  specific 
purpose  of  measuring  Vs  and  Vp  in  dry  sand.  In  this  facility,  a  triaxial 
state  of  stress,  o[  *  o 2  *  o 3  can  be  achieved.  All  tests  performed  to 
date  at  U.T.  have  used  a  local  medium  to  fine,  washed  mortar  sand 
classified  as  SP,  with  effective  grain  size,  Djq  =  0.28  mm  and  a 
uniformity  coefficient  Cu  =  1.7.  The  sand  is  placed  by  the  raining 
technique  and  is  tested  dry.  Values  of  principal  stresses  used  have 
ranged  between  10  psi  and  40  psi,  with  the  stress  ratio,  01/03  =  1  to  4. 
Knox  et  al.  (1982)  used  this  facility  to  study  Vs  and,  similarly  to 
Roesler,  they  propagated  the  shear  waves  along  one  of  the  principal 
stresses  (oa),  and  polarized  parallel  to  another  principal  stress  (ofo). 
They  expressed  Vs  as: 

Vs  =  F  oama  obmb  ocrac  (4) 

where  F  =  constant.  They  found  values  of  ma  =  mb  =  0.09  to  0.12,  and  me 
=  0  to  0.01.  Except  for  some  minor  differences,  these  results  are 
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identical  to  those  of  Roesler,  including  the  independence  of  Vs  on  ac. 
Kopperman  et  al.  (1982)  used  the  same  U.T.  facility  and  sand  to  study 
P-waves  and  concluded  that: 

vp  =  L  oa°*22  (5) 

where  L  =  constant.  The  insensitivity  of  Vp  to  variations  in  the 
stresses  cr^  and  oc  perpendicular  to  wave  propagation  is  illustrated  by 
Fig.  2.  Again,  and  similar  to  Schmertmann's  findings,  these  results 
indicate  that  Vp  is  strongly  dependent  on  direction  of  propagation  when 
the  sand  is  anisotropically  loaded. 

Yu  and  Richart  (1984)  performed  resonant  column  tests  on  three  sands 
subjected  to  a  biaxial  state  of  stress.  Their  results  essentially  agreed 
with  those  of  Roesler  and  Knox;  however,  they  found  some  effect  of  the 
stress  ratio  on  the  results.  They  proposed  for  Gmax  the  expression: 

Gmax  -  CPa°-^9  ov0-26  ah0-25  (1-a  Kn2)  (6) 

where  C=constant,  Pa=atraospheric  pressure,  a=0.15  to  0.23,  with  a  mean 
value  of  0.18,  Kn=(o1/a3~l)/ [(oi/a3)max-l j ,  and  (ai/o2)max  corresponds  to 
shear  failure  of  the  sand.  Except  for  the  factor  1-a  Kn2,  which  is 
usually  between  0.8  and  1,  Eq.  6  is  consistent  with  Eqs.  3  and  4  proposed 
by  Roesler  and  Knox. 

Therefore,  all  of  these  results  clearly  indicate  that  implications 


i)  through  iv)  above,  associated  with  the  currently  used  correlations  for 
Vp  and  Vg  in  sands,  need  to  be  revised  and  upgraded.  In  most  cases  of 
practical  interest,  sands  are  anisotropically  loaded,  and  thus  more  than 
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two  elastic  constants  may  be  needed  to  specify  the  behavior  of  the  soil 
at  very  small  strains.  For  the  typical  biaxial  state  of  stresses  exist¬ 
ing  in  the  field  under  geostatic  conditions,  for  which  all  horizontal 
normal  stresses  are  equal  but  different  from  the  vertical  stress,  the 
sand  will  behave  as  a  cross-anisotropic  elastic  solid  and  5  elastic 
constants  will  be  generally  needed  (Love,  1927,  Sokolnikoff,  1956).  In 
the  more  general  case  of  *  02  *  CT3»  as  ^  may  happen  in  the  soil  under 
a  structure,  the  sand  can  be  described  as  an  orthotropic  elastic  solid, 
with  three  planes  of  elastic  symmetry,  and  a  total  of  9  elastic  constants 
are  needed  (Sokolnikoff,  1956). 

The  previous  discussion  focused  on  the  elastic  properties  of  sand  at 
very  small  strains,  y  *  10 and  especially  on  Vp  and  Vs  measurements. 

If  larger  loads  and  strains  are  applied  to  a  dry  granular  soil,  compres¬ 
sion-wave  type  loading  induces  a  nonlinear  locking  stress-strain  response, 
while  shear-wave  type  loading  induces  a  yielding  response  (see  Fig.  3). 
This  behavior  is  obviously  associated  with  the  particulate  nature  of  the 
soil  (Seed  and  Idriss,  1970;  Hardin  and  Drenvich,  1972).  During  cyclic 
shear  loading  in  sand,  stress-strain  hysteresis  loops  are  generated  such 
as  shown  in  Fig.  4;  these  loops  are  essentially  strain-rate  and  frequency 
independent.  For  small  strains,  y  <  yt  =  the  hysteretic  loop 

repeats  itself  cycle  after  cycle,  and  no  permanent  volumetric  strain  is 
observed,  thus  suggesting  an  essentially  non-destructive  though  nonlinear 
behavior,  controlled  mainly  by  the  response  of  the  contacts  between  the 
grains,  and  with  no  coupling  between  shear  and  volumetric  strains.  At 
shear  strains,  y  >  yt,  although  the  overall  behavior  remains  approximately 
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the  same,  densif ication  occurs,  and  there  is  also  some  increase  in  stiff¬ 
ness,  with  the  shear  stress-strain  curve  and  the  tips  of  the  loop  going 
up  a  little  as  the  number  of  cycles  increases  (see  Fig.  4).  For  the 
strain  range  0.01%  <  y  <  0.1%  to  1%,  the  raonotonic  and  cyclic  behavior  of 
the  sand  is  always  contractive,  that  is,  shear  strains  generate  exclu¬ 
sively  compressive  volumetric  strains,  independently  of  the  density  of 
the  sand.  However,  for  strains  larger  than  about  y  a  1%,  a  mixture  of 
contractive  and  dilative  behavior  is  measured  in  dense  sands,  with 
expansion  of  the  soil  occurring  during  part  of  the  cycle  in  cyclic  shear 
loading  (Youd,  1972).  At  all  strains,  and  for  both  shear  and  compression- 
type  cyclic  loading,  the  stress-strain  response  of  dry  granular  soil  is 
strongly  dependent  on  the  level  of  normal  stresses  acting  on  it. 
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CHAPTER  3 

STRESS-STRAIN  MATHEMATICAL  MODELLING 


A  significant  amount  of  research  has  been  directed  to  obtain 

stress-strain  constitutive  relations  for  cyclic  and  dynamic  loading  of 

soil.  Most  of  these  studies  have  modelled  the  soil  as  an  elastic-plastic 

material,  using  as  a  basis  tool  the  Incremental  Theory  of  Plasticity.  In 

this  type  of  model,  which  is  particularly  appropriate  for  dry  granular 

soil,  the  total  strain  increment  is  equal  to  the  sum  of  elastic  and 

e  p 

plastic  strain  increments,  de  =  de  +  de  ,  with  all  de  being  strain 

ij  ij  ij 

rate  independent  (Drucker  and  Prager,  1952;  Reyes,  1966;  Chen,  1975;  Lade 
and  Duncan,  1975;  Prevost,  1978;  Hardin,  1978).  Based  on  the  Vp  measure¬ 
ments  by  Roesler  previously  described  in  Section  2,  Hardin  (1980)  sug- 

e 

gested  the  following  expressions  for  de  in  dry  granular  soil: 

ij 

dox  doy  do2 

- v  - j. - v  - j 
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F(e) 

pl-n 
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j  e  _  2(1  +  v)F(e)  * 

aYxy  — 


e  e  e 

where  e  =  normal  strain,  y  =  2e  =  engineering  shear  strain,  and  four 
x  xy  xy 

additional  equations  are  obtained  by  permutation  of  subscripts.  In  these 
equations,  crx,  <jy  and  oz  =  normal  stresses;  TXy  =  a  shear  stress;  Pa  = 
atmospheric  pressure;  and  F(e)  =  0.3  +  0.7  e2,  where  e  =  void  ratio. 

Eqs.  7  contain  five  elastic  constants  (Sx,  Sy,  Sz,  SXy  and  v);  based  on 
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Roesler's  experimental  Vs  results,  a  power  of  stress  n  =  0.5  was  pro¬ 
posed. 

A  variety  of  associated  and  nonassociated  flow  roles  have  been 
proposed  for  the  plastic  strain  increment,  of  the  form: 


where  \  =  coefficient  of  proportionality,  and  gCajj)  is  the  plastic 
potential  function,  which  may  or  may  not  coincide  with  the  yield  function 
f(o£j)  at  which  plastic  strains  develop.  Figure  5  shows  the  shapes  of  a 
number  of  plastic  potential  surfaces  proposed  for  soil  by  different 
authors. 

In  the  simplest  type  of  elastic-plastic  model,  there  is  only  one 
yield  (failure)  surface.  For  stresses  below  that  surface,  the  behavior 
is  assumed  to  be  perfectly  elastic.  However,  granular  and  other  soils 
develop  plastic  strains  even  at  the  small  shear  strains  of  interest  to 
this  report.  To  allow  for  this  behavior,  a  wide  variety  of 
strain-hardening  laws  have  been  proposed,  including  families  of  yield 
surfaces  and  specific  strain-hardening  yield  rules.  In  some  of  these 
models,  the  elastic  region  is  completely  eliminated,  thus  allowing  for 
plastic  flow  at  very  low  levels  of  stress  and  strain  (Mroz ,  1967; 

Prevost,  1977).  One  of  the  earliest  developments  included  various  cap 
models,  based  on  the  work  done  at  Cambridge  University  by  Roscoe  and  his 
co-workers  (i.e.,  Roscoe,  1970).  This  includes  the  models  proposed  in 
several  papers  by  DiMaggio  and  Sandler  (1971),  which  have  been  widely 
used  for  dynamic  analyses  of  soil  response  to  explosions.  Several  capped 
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yield  models  are  included  in  Fig.  5. 

An  important  aspect  of  the  development  of  elastic-plastic  models  is 
the  definition  of  the  strain-hardening  law,  which  defines  the  modifica¬ 
tions  of  the  yield  surface(s)  in  the  course  of  the  plastic  flow.  This  is 
especially  critical  for  cyclic  loading,  where  the  type  of  strain-harden¬ 
ing  determines  the  stress-strain  behavior  after  load  reversals.  In  most 
of  the  models  described  above,  which  were  originally  developed  for 
monotonic  loading,  isotropic  strain-hardening  is  assumed  (Hill,  1950), 
with  the  yield  surfaces  expanding  as  the  stresses  increase  (Fig.  6). 

When  isotropic  hardening  is  assumed,  a  large  amount  of  load  reversal  is 
required  for  additional  yielding  to  occur,  in  contradiction  with  the 
observed  behavior  of  experimental  hysteresis  loops  such  as  shown  in  Fig. 
4. 

A  better  alternative  for  earthquake  loading  is  provided  by  the 
kinematic  strain-hardening  law,  sketched  in  Figure  7.  The  kinematic 
model  was  originally  proposed  by  Ishlinsky  (1954)  and  Prager  (1955). 

Iwan  (1967)  proposed  a  reological  representation  for  the  stress-strain 
model,  constituted  by  infinite  elasto-plas tic  elements  placed  in  series 
or  in  parallel.  This  model  is  a  non-f rictional  one,  with  all  the  nested 
yield  surfaces  being  circular  cylinders  in  principal  stress  space.  Mroz 
(1967,  1969)  proposed  a  general  model  for  elastic-plastic  materials, 
composed  also  of  a  field  of  yield  surfaces,  with  a  combination  of 
kinematic  and  isotropic  strain-hardening  laws.  Prevost  (Prevost  and 
Hoeg,  1975;  Prevost,  1977,  1978),  Mroz,  Norris  and  Zienkiewicz  (1978, 
1979)  and  Vicente  and  Dobry  (1983),  have  proposed  to  use  this  model  to 
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predict  static  and  cyclic  behavior  of  soils.  The  model  is  flexible 
enough  to  allow  its  adaptation  to  the  cases  of  drained  and  undrained 
loading,  and  to  incorporate  important  large  strain  cyclic  phenomena  such 
as  densif icatlon,  liquefaction  and  stiffness  degradation.  Anisotropic- 
ally  loaded  soils  are  represented  by  nonsymmetric  nested  surfaces  in 
stress  space.  Under  cyclic  shear  loading,  the  strain-hardening  behavior 
is  basically  kinematic  for  the  reasons  described  above.  A  simultaneous 
isotropic  hardening  (or  softening)  is  allowed  with  the  corresponding 
expansion  or  contraction  of  the  yield  surfaces  as  cyclic  loading 
develops.  This  isotropic  expansion  (contraction)  thus  could  simulate  in 
dry  granular  soil  the  observed  increase  of  stiffness  caused  by  cyclic 
loading  above  a  10~2%. 


CHAPTER  4 


THE  MICROMECHANICAL  APPROACH 

Eqs.  1-2  for  G,,,^  and  similar  expressions  for  other  small  strain 
moduli  in  dry  sands  assume  that  the  controlling  normal  stress  parameter 
is  the  mean  stress  d0:  Gmax  =  f(o0)«  As  discussed  in  Chapter  2,  more 
detailed  measurements  have  revealed  the  elastic  anisotropy  of  a  dry  sand 
subjected  to  different  principal  stresses,  and  they  have  also  shown  that 
the  functional  relationship  between  principal  stresses  a\,  c>2>  a 3,  on  one 
hand,  and  G,,,^  and  other  elastic  constants  on  the  other,  is  not  Gmax  = 
f (oQ)  but  rather  Gmax  »  f(oa,ab>;  similarly,  for  the  constrained  modulus, 

D  =  f(aa).  Derivations  shown  later  in  Section  5.1  for  a  simple  cubic 
regular  array  of  elastic  rough  spheres  match  very  well  with  those  recent 
experimental  findings  in  sands.  This  strongly  suggests  that  a  micro¬ 
mechanical  (particulate  mechanics)  approach  should  be  used  to  analytically 
simulate  and  generalize  the  experimental  observations. 

A  great  number  of  studies  have  been  performed  using  particulate 
models  to  understand  and  model  the  behavior  of  cohesionless  soils  and 
other  granular  materials.  Most  of  these  investigations  have  been  analy¬ 
tical,  but  they  have  also  included  measurements  in  actual  granular  soils, 
as  well  as  in  regular  or  random  arrays  of  spheres  (3-D)  or  disks/rods 
(2-D);  a  number  of  them  have  dealt  with  the  load-deformation  characteris¬ 
tics  at  the  contact  between  two  elastic  bodies  possessing  friction 
(Mindlin's  problem).  State-of-the-art  summaries  have  been  presented  by 
Deresiewicz  (1958,  1973),  Mogami  (1965),  Scott  and  Ko  (1969),  Richart  et 
al.  (1970),  White  (1965),  Harr  (1977),  and  Dobry  and  Grivas  (1978), 
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between  others.  The  proceedings  of  two  US/ Japan  seminars  on  the 
mechanics  of  granular  materials  contain  excellent  papers  on  the  subject 
(Cowin  and  Satake,  1978;  Jenkins  and  Satake,  1983). 

A  number  of  these  studies  have  focused  on  the  probabilistic  aspects 
and  statistical  distributions  of  different  parameters  within  the  soil  or 
granular  medium,  and  their  effect  on  the  mechanical  behavior  of  the 
array.  These  have  included  investigations  of  the  orientations  of  the 
individual  particles,  of  the  spatial  distribution  of  porosity,  and  of  the 
distribution  of  number,  orientation  and  levels  of  force  transmitted  by 
the  contacts,  conducted  between  others  by  Smith  et  al.  (1929),  Dantu 
(1957),  Field  (1963),  Mogami  (1965),  Grivas  and  Harr  (1974),  Oda  (1974), 
Yanagisawa  (1978),  Shahinpoor  and  Shahrpass  (1982),  Nemat-Nasser  and 
Mehrabadi  (1983),  and  Dobry  and  Petrakis  (1984). 

Many  of  those  analytical  investigations,  computer  simulations  and 
observations  have  focused  on  the  stress-strain  behavior  at  very  large 
strains  and  on  the  failure  of  dry  granular  media.  Because  of  this  very 
large  strain  nature  of  the  phenomena,  the  load-deformation  characteris¬ 
tics  of  the  particles'  contacts  have  played  a  minor  or  negligible  role, 
and  the  emphasis  has  been  on  changes  in  the  geometric  arrangement  of  the 
grains  due  to  their  sliding  and  rolling.  In  some  of  these  investigations 
the  compliance  of  the  contacts  has  been  eliminated  altogether  by  assuming 
perfectly  rigid  particles.  Some  important  references  here  are  Rowe 
(1962),  Morgenstern  (1963),  Horne  (1965),  Konishi  (1978),  and  Oda  et  al. 
(1983).  Cundall  and  Strack  (1983)  performed  numerical  experiments  of  2-D 
random  arrays  of  disks  using  an  explicit  finite  difference  procedure 


(see  Fig.  8).  In  these,  the  authors  successfully  simulated  "compression 
triaxial”  loading  to  failure,  and  studied  in  detail  the  spatial  distribu¬ 
tion  of  contact  forces  and  the  distribution  and  relative  contributions  of 


sliding  and  rolling  to  macroscopic  strain,  during  anisotropic  deviatoric 
loading  with  constant  03.  One  of  their  conclusions  for  this  deviatoric 
loading  is  that  the  major  principal  stress  is  transmitted  mainly  by  a 
few  "stiff  chains"  of  particles  having  large  contact  forces,  with  the 
particles  in  between  chains  being  lightly  loaded;  sliding  and  rolling 
occurs  mainly  in  those  lightly  loaded  regions.  Oner  (1984)  worked  with  a 
similar  numerical  scheme  to  predict  the  observed  threshold  strain  yt  at 
which  sliding  and  rolling  starts,  while  Dobry  et  al.  (1982)  used  a 
regular  simple  cubic  array  of  spheres  for  the  same  purpose. 

Of  special  interest  are  the  investigations  which  have  studied  the 
stress-strain  behavior  of  granular  arrays  considering  the  elasticity  of 
the  particles  and  the  corresponding  compliances  at  the  contacts.  Most  of 
these  studies  have  assumed  spherical  grain  shapes  and  elastically  iso¬ 
tropic  grains  characterized  by  three  material  constants  (two  elastic  and 
one  friction  coefficient).  All  of  these  investigations  have  used  the 
normal  and  tangential  compliances  at  the  contact  between  two  elastic 
bodies,  derived  by  Hertz  (1882),  Cattaneo  (1938),  and  Mindlin  and  his 
co-workers  (Mindlin,  1949,  Mindlin  et  al.  1951,  Mindlin  and  Deresiewicz, 
1953).  Figures  9  and  10  show,  respectively,  the  distorsion  of  two 
spheres  subjected  to  normal  (N)  and  tangential  (T)  contact  loads,  and  the 
tangential  load-displacement  curve  for  constant  N.  As  noted  in  the 
summary  reviews  of  the  contact  theory  by  Deresiewicz  (1958,  1973)  and 
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Dobry  and  Grivas  (1978),  Mindlin  and  his  co-workers  developed  the  basic 

theoretical  framework  of  the  contact  problem,  and  solved  it  for  some 

special  force  time  histories;  however,  the  general  problem  of  computing 

+  +  + 

the  displacements  for  a  contact  force  P  =  N  +  T,  whose  magnitude  and 
direction  change  arbitrarily,  remained  unsolved.  Only  very  recently, 

Serldi  and  Dobry  (1984)  provided  a  general  and  practical  solution  to  this 
general  problem,  thus  making  it  possible  the  use  of  direct  stiffness  and 
finite  difference  techniques  to  simulate  the  3-D  response  of  granular 
array  at  small  strains.  A  more  detailed  discussion  of  this  general 
solution  is  presented  in  Section  4.1. 

The  contact  theory  has  been  repeatedly  used  to  predict  the  elastic 
stress-strain  properties  of  granular  arrays  of  spheres.  Several  authors 
calculated  the  influence  of  isotropic  confining  pressure  on  Vp  and  Vs  for 
various  arrays  of  smooth  and  rough  spheres,  and  concluded  that  both 
velocities  increase  proportionally  to  (aQ)^^  (Hara,  1935,  Takahashi  and 
Sato,  1950,  Gassman,  1951,  White  and  Sengbush,  1953,  Brandt,  1955).  Of 
special  interest  here  are  some  detailed  analytical  and  experimental  in¬ 
vestigations  of  regular  arrays  of  equal  spheres.  Deresiewicz  (1958)  lists 
the  five  stable  regular  arrays  included  in  Table  1  and  sketched  in  Fig.  11, 
which  range  from  the  loosest  simple  cubic,  (void  ratio,  e  =  0.91)  to  the 
densest  pyramidal  (also  called  face  centered  cubic  array)  and  tethraedral 
(also  called  hexagonal  close  packed),  both  with  e  =  0.35.  More  complete 
lists  and  descriptions  of  feasible  regular  arrays  have  been  presented  by 
Filep  (1936),  Brown  (1978)  and  Shahinpoor  (1981).  Table  2  reproduces  one 
of  these  lists  containing  31  arrays,  while  Fig.  12  presents  elevation  and 
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plane  views  for  one  of  Che  loosest  arrays  of  Table  2  (Cell  No.  2  with  e  = 
1.94).  Deresiewicz  (1958a)  investigated  in  detail  the  simple  cubic  array 
subjected  to  an  initial  isotropic  loading  followed  by  an  arbitrary  stress 
history.  He  found  that  the  array  is  statically  determinate  in  this  case, 
provided  that  the  stress  field  is  uniform,  with  a  one-to-one  correspon¬ 
dence  between  the  nine  components  of  the  stress  tensor  and  the  nine 
independent  components  of  the  contact  forces.  Whitman  et  al.  (1964) 
studied  a  2-D  version  of  the  simple  cubic  array  subjected  to  triaxial  and 
confined  compression.  Duffy  and  Mindlin  (1957),  Duffy  (1959)  and  Hendron 
(1963)  investigated  the  densest  arrays  of  Fig.  11  and  Table  1,  which  are 
statically  indeterminate,  including  some  measurements  of  compressional 
(rod)  wave  velocity,  Vl,  in  stainless  steel  granular  bars  loaded  iso¬ 
tropically  (see  Fig.  13).  As  shown  in  the  figure,  the  measured  values  of 
are  somewhat  smaller  than  predicted,  with  the  difference  being  greater 
at  small  values  of  cr0,  and  with  this  difference  increasing  for  the  low 
tolerance  balls;  at  high  pressures  the  measured  VL  approach  the  predicted 
one.  As  a  result,  the  observed  a(a0)m,  where  m  >  1/6  =  0.167  pre¬ 
dicted  by  the  theory.  This  difference  is  explained  by  Deresiewicz  (1958), 
by  the  small  differences  in  size  between  the  actual  spheres,  which 
results  in  the  number  of  actual,  load-transmitting  contacts  being  smaller 
than  predicted;  thus,  the  array  is  less  stiff  than  calculated.  When  the 
tolerance  becomes  higher  or  the  pressure  increases,  the  number  of  these 
actual  contacts  also  increases  and  approaches  the  theoretical  value,  and 
thus  the  measured  velocity  also  approaches  the  prediction.  As  discussed 
in  Chapter  2,  values  of  m  a  0.25  >  0.167  have  also  been  measured  for  Vs 
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in  isotropically  loaded  sands,  most  probably  due  to  the  same  reason:  an 


increase  in  the  number  of  actual  contacts  as  a0  increases. 

Several  approaches  have  been  used  to  model  the  effect  of  deviations 
from  regularity  in  arrays  of  equal  or  unequal  spheres.  Smith  et  al. 
(1929)  proposed  considering  a  random  array  as  formed  by  clusters  of  loose 
and  dense  regular  arrays,  each  present  in  such  proportion  as  to  yield  the 
overall  void  ratio  or  porosity;  this  idea  was  generalized  by  Munro  and 
Jowitt  (1974)  and  Brown  (1978),  who  used  the  concept  of  maximum  entropy 
to  find  the  contribution  of  each  regular  array.  Ko  and  Scott  (1967)  used 
a  similar  procedure  to  investigate  the  stress-strain  behavior  under 
isotropic  compression;  in  this  study,  "holey"  models  were  used,  in  which 
some  of  the  spheres  in  both  the  loose  and  the  dense  regular  component 
array  are  slightly  smaller  than  the  other  spheres.  In  this  way  the 
effect  on  the  bulk  modulus  of  the  increased  number  of  contacts  caused  by 
an  increasing  pressure,  was  incorporated  into  the  model.  Perry  and  Brown 
(1981)  studied  the  influence  of  having  different  size  spheres  on  the 
compliance  of  the  array.  Davis  and  Deresiewicz  (1977)  investigated  the 
compressibility  of  a  3-D  random  array  of  smooth  equal  spheres  subjected 
to  isotropic  loading.  Walton  (1987)  derived  a  method  for  computing  the 
effective  elastic  moduli  of  a  random  packing  of  identical  elastic  spheres 
subjected  to  hydrostatic  or  uniaxial  compression.  For  simplicity,  he 
assumed  that  the  spheres  are  either  infinitely  rough  or  perfectly  smooth. 
Serrano  and  Rodrigues-Ort iz  (1973)  suggested  a  method  for  generating 
random  configurations  of  unequal  disks  or  spheres  having  a  prescribed 
grain  size  distribution;  their  work  was  continued  by  the  2-D  numerical 
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simulations  by  Cundall  and  by  Oner,  previously  discussed. 

4.1  General  Solution  to  the  Contact  Problem 

The  problem  of  the  contact  of  two  elastic,  elliptical,  semi-infinite 
bodies  subjected  to  a  normal  force  was  first  studied  by  Hertz  (1882),  who 
also  solved  the  specialized  case  of  the  contact  of  two  elastic  spheres 
subjected  to  a  normal  force,  N.  This  solution  demonstrated  for  the  first 
time  that  the  force-deformation  behavior  at  the  contact  is  nonlinear 
elastic.  Subsequently,  all  work  on  the  same  topic  was  concerned  with  the 
loading  of  bodies  by  normal  forces,  until  Cattaneo  (1938),  Mindlin  (1949), 
and  Mindlin  and  Deresiewicz  (1953)  addressed  the  problem  of  the  contact 
of  two  elastic,  rough  spheres  subjected  to  a  combination  of  normal  and 
tangential  forces,  and  presented  a  number  of  closed  form  solutions  for 
each  of  the  cases  studied.  Recently,  Walton  (1978)  studied  the  problem 
of  the  oblique  contact  of  two  equal  spheres  for  the  case  in  which  both 
the  tangential  normal  forces  are  applied  simultaneously;  this  is  a 
different  case  from  the  work  mentioned  above  and  it  will  not  be  discussed 
here. 

The  general  case  of  two  elastic,  rough  spheres  subjected  to  a  normal 
force  followed  by  a  tangential  force,  (Fig.  9),  solved  by  Cattaneo  (1938) 
and  Mindlin  (1949),  is  a  problem  of  the  linear  theory  of  elasticity; 
however,  since  the  solution  yields  an  infinite  shear  stress  at  the  edge 
of  the  contact  area,  a  slip  needs  to  be  prescribed  at  this  edge,  which 
transforms  the  formulation  into  a  mixed  boundary  value  problem  where  the 
stress  and  the  displacement  are  known  at  the  contact.  This  permanent  set 
induces  a  nonlinear  behavior,  which  is  different  from  that  computed  by 
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Hertz,  since  it  is  accompanied  by  energy  dissipation.  As  demonstrated  by 
Mindlin  and  Deresiewicz  (1953),  due  to  the  slip,  the  force-deformation 
relations  depends  now  on  the  entire  past  history  of  the  loading  as  well 
as  the  instantaneous  rates  of  change  of  the  normal  and  tangential  forces. 
A  typical  force-deformation  behavior  of  two  spheres  subjected  to  a 
constant  normal  force,  N,  and  to  a  monotonically  increasing  tangential 
force,  T,  is  shown  in  Fig.  10,  where  the  nonlinear,  yielding  behavior 
can  be  clearly  observed. 

All  of  the  above  suggest  that  a  phenomenological  plasticity  model 
could  describe  this  nonlinear  behavior.  Such  a  formulation  would  provide 
the  long  awaited  (Detesiewicz ,  1958a)  "general"  solution  to  the  problem 
of  the  contact  of  two  elastic,  rough  spheres  subjected  to  an  arbitrary 
force  history,  which  in  turn,  could  be  used  in  numerical  simulations. 

This  has  been  achieved  recently  through  a  constitutive  law  (Seridi  and 
Uobry,  1984)  for  the  force-deformation  behavior  of  two  identical  elastic, 
rough  spheres  in  contact  under  a  combination  of  arbitrarily  varying 
normal  and  tangential  forces,  which  was  implemented  through  the  program 
CONTACT  (Dobry  et  al.,  1987).  This  model  is  an  elastic-plastic  incremen¬ 
tal  model  with  kinematic  hardening,  the  main  features  of  which  are  similar 
to  those  of  the  plasticity  models  described  in  Chapter  3,  and  are 
described  below: 

Yield  Condition: 

Since  failure  of  the  contact  occurs  when  the  tangential  force, 

2  2  1/2 

T  =  (Tx  +  Ty)  exceeds  fN,  see  Fig.  10,  where  f  is  the  coefficient  of 
friction,  Tx  and  Ty  the  force  compoents  in  the  x  and  y  directions  and  N 
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the  normal  force,  the  failure  surface  is  defined  as  follows: 


T  =  Tx  +  Ty  =  f  N 


This  is  the  equation  for  a  cone  in  the  three-dimensional  force 
space,  N,  Tx,  Ty  (see  Fig.  14). 

Mindlin  and  Deresiewicz  (1953)  demonstrated  that  in  the  loading  case 
where  0  <  dT/dN  <  f,  no  energy  is  dissipated,  while  for  dT/dN  >  f  energy 
is  dissipated.  This  indicates  that  the  yield  surfaces  are  also  cones  of 
semi-angle  <j>  =  tan~^f  parallel  to  the  failure  cone,  the  positions  of 
which  are  determined  by  the  history  of  the  loading,  as  shown  in  Fig.  14. 
Therefore,  the  equation  of  any  yield  surface  is: 

2  2  2  2 

(Tx-x)  +  (Ty-y)  =  f  (N-N£)  (10) 
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where  X  =  x  i  +  y  j  +  N^k  is  the  location  vector  of  the  apex  of  the  conic 

+  A  A  A 

surface,  and  F  =  TX1  +  Tyj  +  N  k  is  the  current  force  point  (Fig.  14). 

Since  the  two  spheres  subjected  to  a  tangential  force  are  constantly 
slipping,  this  implies  that  there  is  an  infinite  number  of  yield  cones. 
The  elastic  region  is  the  inside  of  the  cone  whose  apex  is  at  the  current 
force  point. 

Flow  Rule: 

*  +  + 

For  a  given  force  increment  dF=  dN  -  dT,  the  increment  of  displace¬ 
ment  between  the  centers  of  the  two  spheres  is: 


dD  =  d6x  i  +  d6y  j  +  da  k  =  d6  +  da  k 


The  value  of  dD  is  given  by: 
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dD 


where: 


dN  ft  ,  fdN  *  .  dTn  “  fdN  A  .  dTt  ; 
^_+_n  + - - -  n  +—  t 

1-v 


(12) 


n:  normal  unit  vector  to  the  yield  circle  at  the  current  force 

point, 

A 

t:  tangential  unit  vector  tangent  to  the  yield  circle  at  the 

current  force  point 


dTn  =  dT  .  n 


dTt  =  dT .t  =  (dT2  -  dTn2) 

♦  *  A 

dT  =  dTn  n  +  dTt  t 


1/2 


dTn  +  dTt 


(13) 

(14) 

(15a) 

(15b) 


4G-, 


=  —  is  the  elastic  modulus 

2— v 


,  ,  K  V3 

H  =  Hq  (1  -  -£^)  ,  elastoplastic  modulus  corresponding 


to  the  yield  circle  of  radius  K 


«Ga  K  f  .  _  1  2/3, 

P  2-v  3f N  (1  fN3 


is  the  tangential  elastic  modulus 


(16) 


f:  coefficient  of  friction  of  the  material  of  the  spheres 

a  =*  (BNR)1/3,  with  B  =  -3g3-~V~ 


R  =  radius  of  the  two  spheres 
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Hardening  Rule: 

In  the  3-D  force  space,  the  axis  of  any  yield  surface  translates 
without  rotation  in  such  a  way  that  it  remains  always  parallel  to  the  N- 
axis.  This  is  clearly  a  kinematic  hardening  behavior,  the  mathematical 
expression  of  which  is  given  by: 

+  A  A 

dX  =  dT  -  f  dN  n  +  (K  +  fdN)  dn  (17) 

+ 

where  dX  is  the  translation  in  the  n-plane  of  the  center  of  the  yield 

A 

surface,  K  is  the  radius  of  the  yield  surface  on  the  ir-plane  and  dn  is 
the  change  of  the  direction  in  the  normal  unit  vector  to  the  yield 
surface  on  the  ir-plane  of  the  current  force  point. 

A  sketch  depicting  the  rules  described  above  appears  in  Fig.  14. 

4.1.1  Implementation  of  the  Model  and  Sensitivity  Analysis 

The  above  constitutive  relation  was  implemented  in  the  computer 
program  CONTACT,  as  described  in  detail  by  Seridi  and  Dobry  (1984)  and 
Dobry  et  al.  (1987).  Subsequently,  a  sensitivity  analysis  was  done  in 
order  to  determine  the  range  of  the  values  of  the  force  increment  to  be 
used  in  the  incrementally  linear  analysis.  This  was  done  by  varying  the 
force  increments  in  the  case  in  which  the  normal  force  was  kept  constant 
and  the  tangential  force  increased  monotonically ,  and  by  comparing  the 
resulting  displacements  with  those  obtained  using  the  analytical  solution 
by  Mindlin  and  Deresiewicz  (1953).  It  was  found  that  the  value  of  the 
tangential  displacement  obtained  was  insensitive  to  the  size  of  the  force 
increment,  if  the  ratio  of  the  normal  to  the  tangential  force  increment 
is  less  than  1/50;  it  was  also  found  that  the  error  is  less  than  2 %  if 
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this  ratio  is  about  1/30  (Dobry  et  al.,  1987).  In  all  cases  described  in 
the  next  section,  as  well  as  in  the  numerical  simulations  of  Chapter  7, 
increments  of  similar  size  were  used. 

4.1.2  Verification  of  the  Model 

The  above  constitutive  relation,  implemented  in  the  computer  program 
CONTACT,  was  subsequently  tested  to  verify  that  it  reproduces  accurately 
the  analytical  solutions  obtained  by  Mindlin  and  Deresiewicz  (1953).  Six 
of  the  cases  originally  solved  by  Mindlin  and  Deresiewicz  (1953)  are 
presented  here,  but  this  time  the  force-deformation  behavior  of  the  two 
spheres  in  contact  is  computed  using  the  program  CONTACT  and  the  elastic 
properties  of  quartz  (Es  =  295,182  Kg*/cm2,  vs  =  0.15  and  f  =  0.5,  White, 
1964).  Because  of  the  complexity  of  some  of  the  cases,  only  six  of  the 
most  easily  visualized  cases  will  be  presented  herein.  The  force-defor¬ 
mation  curves  from  the  analytical  results  of  Mindlin  and  Deresiewicz 
(1953)  are  reproduced  separately  in  Figs.  15  and  16,  while  the  correspond¬ 
ing  numerical  results  obtained  with  program  CONTACT  are  shown  in  Figs. 

17  to  22. 

1.  Load-Displacement  Relation  for  Oscillating  Oblique  Force, 

dT/dN  >  f  (Fig.  17): 

The  normal  force,  N,  is  kept  constant  at  1.61  Kg*,  while  the  applied 
tangential  force,  T,  oscillates  between  -0.8  Kg*  and  +0.8  Kg*.  The 
results  are  shown  in  Figure  17,  where  the  data  points  stand  for  the 
analytical  results,  and  the  continuous  curve  for  the  values  obtained 
using  CONTACT. 
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2.  Oscillating  Oblique  Force,  dT/dN  <  f  (Fig.  18): 

The  normal  force  is  increased  to  1.61  Kg*,  and  then  the  tangential 
force  increases  to  0.3  Kg*  while  the  normal  force  decreases  so  that 
dT/dN  =  0.3,  which  is  less  than  f.  Subsequently,  the  shear  force  and  the 
normal  force  oscillate  so  that  dT/dN  remain  constant  at  0.3.  It  can  be 
seen  that  there  is  no  loop  after  one  cycle;  this  was  also  found  by 
Mindlin  and  Deresiewicz  (1953). 


3.  Normal  Force  Increasing,  Tangential  Force  Increasing  (Fig.  19): 

The  normal  force  increases  from  1.61  Kg*  to  2.01  Kg*,  while  the 

tangential  force  increases  from  0.4  to  0.6  Kg*,  and  then  to  1.0  Kg*  (as 
shown  in  the  loading  path  in  the  T-N  space  which  appears  in  Fig.  19). 

The  solid  line  corresponds  to  the  case  in  which  the  tangential  force 
increases  from  0  to  1.0  Kg*  while  the  normal  force  remains  constant.  The 
individual  points  correspond  to  the  force  path  described  above.  This 
corresponds  to  Fig.  7  in  the  Mindlin  and  Deresiewicz  (1953)  paper,  which 
is  reproduced  in  Fig.  16a. 

4.  Normal  Force  Decreasing,  Tangential  Force  Increasing  (Fig.  20): 

The  normal  force  decreases  from  2.01  Kg*  to  1.61  Kg*  while  the 

tangential  force  decreases  from  0.5  Kg*  to  0.3  Kg*  and  then  increases  to 
0.8  Kg*  (as  indicated  by  the  load  path  in  Fig.  20).  The  solid  curve 
corresponds  to  the  force  path  described  above,  and  the  individual  points 
to  the  case  in  which  the  tangential  force  increases  with  the  normal  force 
remaining  constant  (this  corresponds  to  Fig.  16b,  which  is  Fig.  10  in 
Mindlin  and  Deresiewicz  (1953)). 
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5.  Normal  Force  Increasing,  Tangential  Force  Decreasing  (Fig.  21): 

The  normal  force  increases  from  1.61  Kg*  to  2.01  Kg*,  while  the 

tangential  force  is  reduced  from  0.4  to  0.2  Kg*  (as  indicated  in  Fig.  21). 
Again,  the  solid  lines  correspond  to  the  stress  path  described  above,  and 
the  individual  points  to  the  force  path  in  which  the  tangential  force 
increases  from  0  to  1.0  Kg*  while  the  normal  force  remains  constant  at 
1.61  Kg*  (this  corresponds  to  Fig.  16c,  which  is  Fig.  12  in  Mindlin  and 
Deresiewicz  (1953)). 

6.  Normal  Force  Decreasing,  Tangential  Force  Decreasing  (Fig.  22): 

This  is  a  complicated  loading  path  (shown  in  Fig.  22).  The  points, 

once  more,  correspond  to  the  force  path  with  the  constant  normal  force. 
This  reproduces  Fig.  14  of  Mindlin  and  Deresiewicz,  which  has  been 
reproduced  in  Fig.  16d. 

It  can  be  seen  that  the  computer  program  CONTACT,  which  is  based  on 
the  constitutive  law  derived  by  Seridi  and  Dobry  (1984),  reproduces 
accurately  all  the  solutions  obtained  by  Mindlin  and  Deresiewicz.  There¬ 
fore,  it  can  be  assumed  with  confidence  that  the  program  functions 
properly.  On  this  basis  it  will  be  used  in  the  numerical  simulations 
of  Chapter  7. 


4.2  Studies  in  Crystals 

The  concept  of  Micromechanics,  besides  its  obvious  application  to 


the  behavior  of  granular  media,  has  been  applied  over  the  years  to  pre¬ 
dict  the  response  of  polycrystalline  aggregates  from  the  behavior  of 
single  crystals.  It  has  been  demonstrated  that  this  micromechanical 


modeling  can  adequately  interpret  most  of  the  phenomena  associated  with 
the  elastic  and  plastic  macroscopic  behavior  of  polycrystalline  media 
such  as  metals. 

For  years,  researchers  have  been  simulating  the  elastoplastic 
behavior  of  polycrystalline  aggregates  through  a  variety  of  analytical, 
semianalytical  and  numerical  techniques  (Taylor,  1938,  Hershey,  1954, 
Bishop  and  Hill,  1951,  1951a,  Budiansky  et  al»,  1960,  Budiansky  and  Wu, 
1962,  Lin,  1964,  Lin  and  Ito,  1966,  1967,  Hill,  1967,  Canova  et  al., 

1984,  1985).  One  of  the  approaches  most  commonly  used  is  to  model  the 
polycrystal  as  an  assemblage  of  equal,  isotropic  crystals,  with  N  slip 
systems,  which  are  randomly  oriented  (Fig.  23).  This  results  in  an 
isotropic  polycrystal  when  the  distribution  of  the  orientation  of  the 
crystals  is  statistically  uniform;  moreover  this  distributes  randomly  in 
all  directions  the  orientations  of  the  slip  planes  of  the  monocrystals . 
The  yield  surfaces  of  the  resulting  macroscopic  medium  are  determined  by 
these  slip  planes.  Plastic  strain  in  the  aggregate  starts  when  the  first 
individual  crystal  slips;  after  sliding  has  occurred  in  a  number  of  these 
monocrystals,  the  yield  surface(s)  change  (harden)  accordingly.  As  the 
aggregate  is  loaded  farther  beyond  the  elastic  range,  more  and  more 
crystals  slip,  and  eventually  the  failure  surface  is  reached  (Lin  and 
Ito,  1966). 

In  the  early  stages  of  that  research,  Taylor  (1938),  by  considering 
an  aggregate  of  randomly  oriented  rigid-plastic  for  crystals  under 
uniaxial  tension,  computed  the  amounts  of  slip  in  each  crystal  by 
assuming  that  the  strain  experienced  by  each  monocrystal  was  the  same  as 


the  applied,  macroscopic,  uniform  strain.  Then,  by  equating  the  sum  of 
the  work  in  the  monocrystals  to  the  work  done  on  the  aggregate,  he  was 
able  to  obtain  a  tensile  stress-strain  curve,  and  a  value  for  the  yield 
stress  of  the  aggregate  which  later  was  proved  to  be  the  upper  bound  of 
the  actual  yield  stress  (Lin,  1964).  Taylor  also  showed  that  the  slip  in 
the  monocrystal  depends  on  the  resolved  shear  stress  (the  shear  stress 
along  the  slip  direction  and  on  the  slip  plane),  and  that  this  slip  is 
independent  of  the  normal  stress  acting  on  that  plane.  Bishop  and  Hill 
(1951)  derived  a  relationship  between  the  stress  and  the  plastic  strain 
in  a  polycrystalline  aggregate  of  fee  crystals,  and  obtained  a  more  re¬ 
fined  value  of  the  limiting  value  of  the  yield  stress  which  they  found  to 
be  1.536ay,  where  oy  is  the  yield  stress  of  the  single  crystal.  The  same 
valu»  fas  computed  by  Budiansky  et  al.  (1960),  who  abandoned  the  assump¬ 
tion  of  strain  uniformity  throughout  the  polycrystal  and  computed  instead 
the  stress  field  inside  each  crystal.  Using  general  results  of  the 
Eshelby  (1957)  solution  for  the  elastic  field  inside  an  elliptical,  elas¬ 
tic  inclusion  contained  within  an  elastic  medium,  they  solved  the  problem 
of  an  isolated,  elastic-ideally  plastic,  spherical  grain  contained  within 
an  infinite,  stressed  elastic  matrix,  and  they  proved  that  the  stress  and 
strain  fields  inside  the  grain,  while  orientation  dependent,  were  homo¬ 
geneous.  Later,  Budiansky  and  Wu  (1962)  considered  an  "average"  inter¬ 
action  between  the  different  phases,  that  is  the  effect  of  the  stress 
release  due  to  a  slip  in  one  group  of  crystals  on  the  corresponding 
average  stress  increase  in  another,  and  they  calculated  numerically 
macroscopic  stress-strain  curves  in  simple  tension  and  pure  shear. 


The  behavior  of  a  random  packing  of  identical  spheres  is  similar  in 
some  respects  to  the  polycrystalline  medium  just  described,  as  it  is  well 


known  that  random  assemblages  of  equal  spheres  tend  to  crystalize  (see 
Chapter  6).  Since  the  individual  grains  in  a  uniform,  granular,  rounded 
grained  soil  can  be  modelled  by  spheres,  it  is  possible  to  consider,  in 
first  approximation,  that  the  individual  grain  packings  within  the  sand 
behave  like  randomly  oriented  crystals.  However,  a  main  difference  is 
that  the  properties  of  these  grain  packings  are  now  pressure  dependent 
and  the  amount  of  slip  in  each  of  these  packings,  in  contrast  to  the 
polycrystalline  aggregates,  depends  now  on  the  normal  pressure  acting  on 
the  slip  plane.  For  example,  it  can  be  assumed  that  a  simple  cubic  array 
of  equal  spheres  (sand  grains)  is  a  pressure  dependent  monocrystal  with 
three  slip  planes,  with  each  slip  plane  containing  two  slip  directions  90 
degrees  apart.  Also,  at  large  strains  beyond  the  first  slip  (correspond¬ 
ing  to  the  threshold  strain  of  the  soil),  the  sand  may  experience 
dilation  under  shear  which  is  not  necessarily  present  in  polycrystalline 
aggregates. 


I 
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CHAPTER  5 

DIFFERENTIAL  STRESS-STRAIN  RELATIONS  FOR  REGULAR  ARRAYS  OF  SPHERES 

The  general  solution  of  the  problem  of  the  contact  between  two 
spheres  can  be  used  to  derive  incremental  stress-strain  relationships  for 
regular  arrays  of  spheres.  These  stress-strain  relationships  are 
discussed  in  this  Chapter,  with  particular  emphasis  on  the  behavior  under 
isotropic  loading  followed  by  very  small  but  arbitrary  stress  and  strain 
increments,  and  for  the  following  regular  arrays: 

i)  simple  cubic  array  (sc,  see  Figs.  11a  and  24),  discussed  in 
Section  5.1 

ii)  body  centered  cubic  array  (bcc,  see  Fig.  25);  discussed  in 
Section  5.2 

iii)  face  centered  cubic  or  pyramidal  array  (fee,  see  Figs,  lid  and 
26);  discussed  in  Section  5.3 

iv)  cubical  tetrahedral  (ct,  see  Fig.  lib)  and  tetragonal-sphenoidal 
arrays  (ts,  see  Fig.  11c);  discussed  in  Section  5.4. 

Most  of  these  incremental  stress-strain  relations  were  taken  from 
Deresiewicz  (1958a),  Duffy  and  Mindlin  (1957)  and  Maklhouf  and  Stewart 
(1967).  However,  some  are  new;  in  particular,  the  body  centered  cubic 
array  is  discussed  here  for  the  first  time. 

In  addition  to  the  five  regular  arrays  listed  above,  a  sixth  regular 
array  is  the  hexagonal  closed  packed  or  tetrahedral  array  (hep,  see  Fig. 
lie,  see  also  Deresiewicz,  1958).  Table  1  lists  the  most  important 
parameters  of  all  six  arrays.  A  comparison  of  the  behavior  of  the  sc, 


bcc  and  fee  cubic  arrays  is  presented  in  Section  5.5. 

5.1  Simple  Cubic  Array  (sc) 

The  simple  cubic  array  sketched  in  Figs.  11a  and  24  is  the  simplest 
of  all  regular  arrays  of  equal  spheres.  One  sphere  of  radius  R  represents 
the  whole  array,  and  for  a  uniform  stress  field  this  is  a  statically 
determinate  system  with  a  one-to-one  correspondence  between  the  array's 
stresses  and  the  contact  forces  (Deresiewicz ,  1958a).  If  the  normal 
stresses  parallel  to  the  axes  of  the  array  are  (i  =  1,2,3,  see  Fig. 
24),  the  normal  contact  forces  are:  =  4R^  If  the  shear 

stresses  parallel  to  the  axes  of  the  array  are  oij  (i,j  =  1,2,3  and  i  t 
j),  the  corresponding  tangential  contact  forces  and  T^j  =  4R^  j .  These 
relations  occur  due  to  equilibrium  and  are  independent  of  the  previous 
history  of  stresses.  Therefore,  they  are  also  valid  for  any  stress  and 
force  increments  at  any  stage  during  the  loading,  provided  that  no  gross 
sliding  of  the  contact  has  taken  place  dN^  =  4R^  aii,  dTjj  =  4R^  do-jj. 
Figure  24  illustrates  the  case  in  which  an  anisotropic  state  of  stresses 
is  applied  first,  with  all  a^j  =  0,  followed  by  small  arbitrary 
increments  do^  and  do^ j . 

A  similar  set  of  simple  relations  is  valid  in  this  case  between  the 
array  strains  and  the  displacements  between  spheres;  these  relations  are 
obtained  for  a  uniform  strain  field  based  on  simple  geometric  considera¬ 
tions.  If  cii  *  normal  relative  displacement  of  centers  of  the  two 

Indicial  tensor  notation  is  not  used  here,  that  is,  cr^i  does  not 
imply  a  sum  of  several  terms. 
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adjacent  spheres  separated  by  contact  i,  and  5ij.  6ilc  =  tangential  rela¬ 
tive  displacements  between  the  same  two  centers,  then:  En  =  oti/(2R); 
and  Yij  =  2eij  =  (6ij  +  6ji)/(2R),  where  en  =  normal  strain,  and  Yij  = 
engineering  shear  strain  of  the  array.  Again,  all  these  relations  are 
independent  of  the  prior  history  of  strains  and  are  valid  for  incremental 
displacements  and  strains. 

The  elastic  stiffnesses  corresponding  to  small  stress  and  strain 
increments  applied  to  the  array  subsequent  to  an  isotropic  stress  state, 
CH  *  a22  =  a23  =  aot  0ij  =  o  are: 
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Notice  that  the  stiffness  matrix  is  diagonal,  and  therefore  the 
Poisson's  ratio  of  the  array,  v  =  | de22/de 1 l I  =  for  "triaxial"  loading 
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corresponding  Co  increasing  on  and  constant  022  ~  a33  ~  °om  array 

has  a  v  =  0  quite  independently  of  the  values  of  o0  and  vs  (see  Fig. 
29). 


Note  that,  for  a  given  a0,  Eqs .  18-19  describe  a  linearly  elastic 
anisotropic  medium.  The  necessary  and  sufficient  conditions  for  this 
medium  to  be  isotropic  are  of  the  type  Sjjji  -  S{{22  =  ^1212*  These  are 
satisfied  only  for  a  Poisson's  ratio  of  the  spheres,  vs  =  0.  This 
assumption  results  in  an  Isotropic  medium  with  v  =  0. 

For  the  case  of  an  anisotropic  state  of  stress,  on  *  022  *  033,  o^j 
=  0,  (i  +  j),  the  coefficients  become: 
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where  Gg  and  vg  are  the  shear  modulus  and  Poisson's  Ratio  of  the  material 
of  the  spheres. 

The  array  locks  under  "triaxial”  conditions  (o^  increasing  with 
022  “  033  =  0o  =  constant),  while  it  fails  in  pure  shear  (012  increasing 
with  all  other  do  =  0).  Because  the  stiffness  matrix,  [S] ,  in  Eq .  18  is 


diagonal,  the  corresponding  diagonal  compliance  matrix  is  [CJ  =  [S]~^, 
with  each  compliance  being  the  reciprocal  of  the  corresponding  diagonal 
stiffness  term.  That  is,  =  1/S^j^  and  Cjjjj  a  1/Sijij,  and 

{de}  =  [C]  {da}. 

In  this  simple  cubic  array,  and  again  for  the  case  of  isotropic 
loading,  =  ajj  =  0^=  aQ;  for  P-  and  S-waves  propagating  along  prin¬ 

cipal  axis  i  of  the  array  (and,  for  the  S-wave,  polarized  parallel  to 
principal  axis  j),  the  wave  propagation  velocities  Vp  and  Vs  are  propor¬ 
tional  to  (a0)*/6,  and  thus,  the  corresponding  constrained  and  shear 
,  2 

moduli,  D  =  pVp^  =  do^/deii  and,  G^x  =  pVs  =  daij/dvij,  are  both  pro¬ 
portional  to  (ao)0*33.  As  discussed  in  Section  4,  a  similar  dependence 
of  modulus  on  is  also  predicted  for  other  regular  arrays,  while 

laboratory  measurements  on  regular  arrays  and  soils  indicate  that  Gmax 
a(o0)0,5. 

For  this  same  case  of  isotropic  loading  and  for  a  cubic  array  of 
quartz  spheres,  and  if  the  array  is  loaded  in  pure  shear,  dajj  =  doji, 
then  a  threshold  shear  strain,  =  4.5  x  lO-^.  (a0)^/3  ts  predicted 
where  Yt  is  in  inch/in  and  aQ  is  in  psi  (see  Appendix  B,  Eq .  B5).  This 
expression  was  obtained  using  the  properties  of  quartz  listed  in  Table  3. 
At  y  =  Yt»  gross  sliding  occurs  at  contacts  i  and  j,  and  there  is  a 
tendency  for  a  change  in  the  geometric  arrangement  of  the  spheres. 

This  predicted  relation  between  yt  and  aQ  for  a  sc  array  of  quartz 
spheres  is  plotted  in  Fig.  22  while  Fig.  28  shows  the  detailed  shear 
stress-strain  plots  up  to  the  threshold  (failure).  In  the  range  of  most 
practical  interest  for  soils,  5  psi  <  o0  <  15  psi,  the  expression  gives 
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yt  =  10 ~^%  =  10~2%,  which  agrees  very  well  with  the  measured  yt  sands 
as  discussed  in  Section  2.  Expressions  for  the  secant  modulus  reduction, 
G/Gmax»  and  for  the  damping  ratio  of  the  array,  versus  strain  increment  y 
=  were  also  obtained  for  a  cubic  array  of  quartz  spheres  (Dobry  et 

al.,  1982),  and  were  compared  with  actual  measurements  in  sands,  with  good 
agreement.  The  corresponding  comparison  for  G/Graax  versus  y  is  reproduced 
in  Fig.  29  for  an  assumed  =  1.5  x  10 ~2%. 

The  more  general  case  of  anisotropic  loading  of  the  cubic  array, 
with  ajj  *  022  *  a33»  very  interesting,  as  this  model  crudely  simu¬ 
lates  the  laboratory  measurements  of  Vp  and  Vs  on  anisotropically  loaded 
sands  discussed  in  Chapter  2.  For  a  P-wave  propagating  parallel  to  the 
i-axis  of  the  array,  the  predicted  expressions  for  D  =  dan/den  and  Vp 
are: 


2  2/3  1/3 

D  -  [(3)1/3/2][Es/(1-vs)  Koh) 

Vp  =  (D/p)l/2 


(22) 


where  Eg,  vs  =  elastic  constants  of  the  spheres.  Therefore,  both  D  and 
Vp  are  functions  only  of  the  normal  stress  an  in  the  direction  of 
propagation,  and  do  not  depend  on  the  other  two  array  stresses  ojj  and 
°kk* 


For  an  S-wave  propagating  parallel  to  the  i-axis  of  the  anisotropi¬ 
cally  loaded  array  and  with  motions  polarized  parallel  to  the  j-axis  of 
the  array,  the  corresponding  expressions  for  G^*  =  dajj/dy^j  and  Vg  are: 


Gmax  and  Vs  are  functions  only  of  the  normal  stresses  in  the  direction  of 
propagation  (a^)  and  polarization  (ojj),  and  do  not  depend  on  the  third, 
out  of  plane  array  stress  okk*  Furthermore,  as  Eqs.  23  are  symmetric 
with  respect  to  an  and  ajj,  the  values  of  Vs  and  G,,,^  do  not  change  if 
the  directions  of  propagation  and  polarization  are  interchanged. 

These  conclusions  for  the  simple  cubic  array,  that  Vp  depends  only 
on  an,  and  Vs  depends  only  on  a  a  and  ojj,  are  identical  to  the  experi¬ 
mental  findings  of  Roesler  (1979),  Knox  et  al.  (1932),  Kopperman  et  al. 
(1982),  and  Yu  and  Richart  (1984)  on  anisotropically  loaded  sands, 
previously  discussed  in  Section  2.  The  symmetry  of  on  and  ojj  in  Eq.  23 

is  also  present  as  a  symmetry  of  oa  and  Ob  in  empirical  Equation  4, 

obtained  by  Knox  et  al.  (1982)  from  their  sand  measurements. 

Of  course,  Eqs.  22  and  23  cannot  be  used  directly  for  quantitative 
predictions  in  sands,  as  they  give  D  a  and  Gmax  a  o0 for  the 

isotropic  case,  while  D  and  C^jax  a(o0)^^  in  actual  sands.  It  is  inter¬ 
esting  to  modify  Eq.  23  to  make  it  consistent  with  this  empirical  fact, 

by  replacing  ojj^^  by  and  then  to  compare 

measurements  and  predictions.  The  new  equation  is  -’oj  j®*  V 

(oiiO-S+ojjO-S),  where  M  =  constant.  It  is  useful  to  specialize  this 
expression  for  the  biaxial  case,  o^£  =  oil  =  ov,  ojj  =  o^k  =  °33  =  ah*  t0 
be  able  to  compare  it  with  the  empirical  Eq.  6  obtained  by  Yu  and  Richart 
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for  sand.  If  K  =  ^ii/ojj  =  oil^°33>  the  new  equation  becomes:  G^x  = 
Nav^* 25[  2j(0. 25/ (i+^O.  5)  ]  f  where  N  =  0.5M.  It  is  convenient  to  de¬ 
fine  the  normalized  parameter  G  =  (Nov® •  25^0. 25)  f  where  G  =  1  for 

K  =  1.  The  theoretical  expression  G  =  2kP*25/(1+kO*5)  has  been  plotted 
in  Fig.  30.  The  corresponding  empirical  expression  G  =  1-0.18 
obtained  from  Eq.  6,  has  also  been  superimposed  on  Fig.  30  for  typical 
values  Kgjax  =  3  and  4.  The  trends  of  the  predicted  and  measured  curves 
are  the  same  in  Fig.  30,  with  the  laboratory  results  showing  a  somewhat 
faster  decrease  In  G  as  K,,  increases. 

The  fact  that  the  crude  particulate  model  used  here  is  capable  of 
predicting  the  lack  of  influence  of  the  two  stresses  perpendicular  to 
propagation  on  Vp  (Eq.  22),  and  of  the  out-of-plane  stress  on  Vg  and  G^ax 
(Eq.  23),  as  well  as  the  general  trend  of  the  relationship  between  Gmax 
and  the  in-plane  stresses  (Fig.  30),  is  extremely  encouraging.  The  main 
advantage  of  the  cubic  array  used  here  is  its  simplicity,  but  of  course 
this  model  is  still  far  from  representing  real  sand.  One  deficiency 
(which  it  shares  with  other  regular  arrays),  is  that  in  the  general  case 
the  array  itself  is  inherently  anisotropic  even  when  isotropically  loaded 
(crystal-type  behavior);  that  is,  Graax  and  other  "elastic"  stress-strain 
parameters  are  somewhat  different  for  shear  stresses  corresponding  to 
axes  which  are  different  from  the  structural  axes  (axes  of  symmetry  of 
the  array  1,  2,  3  selected  in  Fig.  24).  However,  when  the  material  of 
the  spheres  has  vs  =  0,  the  array  is  isotropic  when  isotropically  loaded 
at  very  small  increments.  Also,  the  array  locks  when  a  "triaxial”  test 
is  conducted  on  it  along  any  of  its  structural  axes,  instead  of  yielding 
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and  eventually  falling  in  shear  as  it  happens  with  actual  granular 
materials. 


5.2  Body  Centered  Cubic  Array  (bcc) 

The  body  centered  cubic  array,  sketched  in  Fig.  25,  was  the  next 
regular  array  studied.  It  is  also  represented  by  one  sphere,  and  the 
relations  between  stress  and  contact  forces  can  be  easily  determined  for 
one  uniform  stress  field  of  interest:  isotropic  loading  followed  by 
small  stress  increments.  The  coordination  number  (number  of  contacts/ 
sphere)  is  now  eight  instead  of  six  for  the  simple  cubic  array,  and, 
thus,  the  computations  are  somewhat  more  involved.  A  procedure  analogous 
to  that  used  for  analyzing  the  simple  cubic  array  can  be  followed,  except 
that  it  is  now  easier  to  work,  directly  with  the  compliance  matrix,  Cijkfc* 
[C]  =  [S]~l,  instead  of  the  stiffness  matrix  [S]  used  before  in  Eq.  18. 

For  the  case  of  isotropic  loading,  followed  by  small  stress 
increments,  [C]  has  the  following  form: 
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where 
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Notice  that,  unlike  Eq.  18  for  the  sc  array,  the  compliance  matrix  [C]  in 
Eq.  24  is  not  diagonal.  However,  it  becomes  diagonal  and  it  corresponds 
to  an  isotropic  elastic  medium  with  v=  0,  if  vs  =  0,  similarly  to  that 
found  for  the  sc  in  Section  5.1.  For  an  initial  cross-anisotropic  or 
biaxial  loading,  022  =  oQ  +  oa  and  on  =  033  =  a0,  followed  by  arbitrary, 
small  stress  increments*,  the  form  of  [C]  is  still  that  of  Eq.  24,  and 


the  compliances  in  Eq.  24  are: 
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See  footnote.  Appendix  A. 3. 
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As  it  can  be  seen  in  the  above  equations  25-30,  Vs  and  Vp  are  again  pro¬ 
portional  to  (<j0)1/6  for  the  body  centered  cubic  array,  since  the  cor¬ 
responding  moduli  are  proportional  to  (o0)^3;  this  is  a  characteristic 
common  to  all  regular  cubic  arrays  (Duffy  and  Mindlin,  1957,  Duffy,  1959, 
Makhlouf  and  Stewart,  1967).  However,  again  it  is  possible  to  modify  the 
exponents  empirically,  and  (oq)^3  can  be  replaced  by  (o0)^'2  when 
measurements  and  predictions  are  compared. 

Threshold  strain  calculations  were  performed  for  this  body-centered 
array  in  the  case  of  triaxial  loading,  starting  from  an  isotropic  pres¬ 
sure  o0,  for  which  the  array  yields  and  fails  (while  the  array  tends  to 
lock  under  pure  shear  loading).  Again,  the  threshold  strain,  yt ,  obtained 
for  the  array  is  a  function  of  the  confining  pressure,  (o0)2/3,  and  for 
an  array  of  quartz  spheres  (using  the  properties  for  quartz  in  Table  3), 

Yt  “  3.44  x  1 0-3  gC)2/3>  with  oQ  in  psi  and  yt  in  inches/inch.  This  gives 
slightly  lower  values  of  yt  than  those  obtained  from  the  simple  cubic  ar¬ 
ray  of  quartz  spheres  in  Section  5.1,  yt  =  4.53  x  ID-3  gQ2/3.  The  plots 
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for  these  two  expressions  of  yc  are  compared  in  Fig.  27,  while  Fig.  31a) 
presents  detailed  axial  stress-strain  curves  up  to  the  threshold  (failure), 
For  the  usual  range  of  values  of  confining  pressure  for  soils,  both  of 
them  agree  well  with  yt  =  10 experimentally  observed  in  sands. 

The  body  centered  cubic  array  also  has  some  deficiencies  when  com¬ 
pared  to  the  behavior  of  actual  sands.  First,  it  remains  isotropic  even 
when  loaded  under  anisotropic  loads,  as  shown  in  Appendix  A.  Second,  in 
this  anisotropic  loading  case  the  wave  propagation  velocities  are  not 
proportional  to  the  product  of  the  principal  stresses  in  the  directions 
of  propagation  and  polarization  (as  measured  in  sands),  but  rather  they 
are  proportional  to  the  mean  stress,  as  can  be  verified  from  Eqs.  28-30. 
Finally,  this  array  locks  under  pure  shear  loading,  and  in  fact  it  locks 
under  a  number  of  different  shear  loading  paths  depending  upon  their 
orientations  and  initial  stress  state.  However,  the  bcc  array  adds  to 
our  understanding  of  the  general  problem,  as  it  is  a  medium  dense  (e  = 
0.47)  array,  located  within  the  range  between  the  densest,  face  centered 
cubic  array  (e  =  0.35)  and  the  loosest,  simple  cubic  array  (e  =  0.91). 


5.3  Face  Centered  Cubic  Array  (fee) 

The  Face  Centered  Cubic  Array  sketched  in  Fig.  26  is  one  of  the  two 
densest  arrays,  and  it  has  been  investigated  by  several  researchers 
(Duffy  and  Mindlin  1957,  Ko  and  Scott  1969,  Hendron  1963).  The  differen¬ 
tial  stress-strain  relationship  for  this  medium  was  derived  by  Duffy  and 
Mindlin  (1957).  The  array  has  12  contacts  per  sphere,  and  unfortunately 
it  is  statically  indeterminate  for  most  loading  situations;  as  a  conse¬ 
quence,  closed  form  solutions  are  available  only  for  the  c  ase  of 


isotropic  confining  pressure.  In  the  case  of  biaxial  loading  a  qualita¬ 
tive  solution  does  exist,  but  the  compliances  at  the  contacts  have  not 
yet  been  evaluated.  Computation  of  these  compliances  is  a  formidable 
task  due  to  the  indeterminacy  of  the  problem  and  the  variation  of  the 
forces  from  contact  to  contact. 


The  incremental  constitutive  law  under  isotropic  confining  pressure 
was  used  by  Duffy  and  Mindlin  to  compare  the  theoretical  and  experimental 
rod  wave  velocities  through  a  bar  composed  of  face  centered  cubic  arrays 
of  spheres  (Fig.  13). 

For  the  case  of  isotropic  loading  followed  by  small  stress 
increments,  don,  dojj ,  stiffness  matrix  has  the  form  shown  below: 
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Similar  to  that  discussed  previously  for  Eq.  24  and  the  bcc  array,  the 
stiffness  matrix  [S]  in  Eq.  31  becomes  isotropic  and  diagonal  only  if  vs 
=  0.  In  that  case,  all  diagonal  terms  are  identical  and  equal  to: 


3<^  a0 
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(35) 


and  the  Poisson's  ratio  of  the  array  is  v  =  0. 

If  the  face  centered  cubic  array  is  consolidated  under  a  transversely 
isotropic  state  of  stress,  with  oji  *  aQ  +  oa,  022  =  °33  =  °o>  and  °12  = 
013  *  °23  “  0,  the  situation  becomes  more  complicated,  since  in  the  case 
of  a  non  isotropic  loading  the  forces  vary  from  contact  to  contact  and 
each  compliance  is  different.  To  obtain  a  stress-strain  relation  for  an 
anisotropic  loading  path,  the  derivation  must  be  performed  anew,  distin¬ 
guishing  between  contacts  with  different  loading  histories. 

The  differential  stress  strain  law  for  this  case  appears  in  great 
detail  in  the  original  paper  by  Duffy  and  Mindlin  (1957);  it  is  identical 
in  form  to  Eq.  (31),  except  that  now  the  expressions  for  are  not 

known.  Thurston  (1958)  extended  the  results  of  Duffy  and  Mindlin  to  a 
set  of  18  equations  and  18  unknowns. 

The  fee  array  was  subjected  analytically  to  the  conditions  of  a 
triaxial  test  by  Brauns  (1968)  and  Brauns  and  Leussink  (1970),  who 
derived  theoretical  expressions  between  stress  and  strain  at  finite 
levels  for  an  array  of  glass  spheres  (Fig.  31b).  These  expressions  were 
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later  compared  to  experimental  data  obtained  in  triaxial  tests  on  regular 
fee  packings  of  glass  and  steel  spheres  (see  Appendix  B3). 

Thurston  and  Deresiewicz  (1959)  derived  expressions  for  the  uniaxial 
compression  of  an  fee  array  when  applied  concurrently  with  a  related 
isotropic  pressure.  Again,  the  theoretical  results  were  compared  with 
experimental  results  obtained  through  compression  of  bars  of  steel 
bearing  spheres  arranged  in  fee  array,  with  good  agreement. 


(37) 
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As  we  can  see  from  the  above  equations,  the  Cubical  Tetrahedral  array 
differs  from  the  simple  cubic,  the  body  centered  cubic  and  the  force 
centered  cubic  arrays  in  that  it  does  not  exhibit  cubic  anisotropy,  but 
rather  transverse  or  hexagonal  isotropy,  as  is  also  the  case  of  the  Hexa¬ 
gonal  Close  Packed  array  (Duffy,  1959)  and  of  the  Tetragonal-Sphenoidal 
Array. 

The  constitutive  law  for  the  Tetragonal-Sphenoidal  array  appears 
also  in  Makhlouf  and  Stewart  (1967).  However,  not  enough  detail  is 
provided  in  this  reference  for  a  full  underst u  ! f ng  of  the  results,  which 
are  quite  complex,  as  the  representative  unit  prism  is  not  symmetric. 
Unfortunately,  the  original  reference  (Makhlouf,  1963)  could  not  be  found 
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by  the  authors*,  thus  preventing  a  better  understanding  of  this  array. 

5.5  Comparison  of  Different  Cubic  Arrays 
The  analytical  results  on  the  three  regular  cubic  arrays  discussed 
in  Sections  5. 1-5. 2-5. 3:  simple  cubic,  body  centered  cubic  and  face 
centered  cubic,  were  compared  as  part  of  the  current  research.  This  was 
done  to  gain  further  insight  into  the  behavior  of  granular  media,  and  as 
a  necessary  intermediate  step  toward  the  investigation  of  more  elaborated 
and  realistic  particulate  models. 

All  these  arrays  generally  exhibit  cubic  anisotropy  (crystal-type 
behavior)  under  an  isotropic  confining  pressure  o0.  In  the  three  arrays, 
it  was  found  that  the  necessary  and  sufficient  condition  for  the  array  to 
become  isotropic  under  oQ  is  for  the  Poisson's  ratio  of  the  spheres,  vs, 
to  be  equal  to  zero.  If  vs  *  0,  the  incremental  stiffness  (and  compli¬ 
ance)  matrix  for  the  three  arrays  is  diagonal.  Although  the  Poisson's 
Ratio  of  quartz  vs  »  0.15,  (see  Table  3),  is  certainly  different  from 


The  differential  constitutive  laws  for  the  cubical  tetrahedral  and  the 
tetragonal  sphenoidal  arrays  are  either  not  applicable  to  our  research  or 
they  are  erroneous.  As  one  can  see  from  Eqs.  36-40,  contrary  to  general 
belief  (Duffy  1959),  the  cubical  tetrahedral  array  does  not  become  iso¬ 
tropic  under  isotropic  loading.  This  is  a  serious  deficiency  vis-a-vis 
our  research,  as  sands  are  uniform  and  more  or  less  isotropic  under 
isotropic  load.  Consequently,  the  cubical  tetrahedral  array  will  not  be 
used  here.  As  for  the  tetrogonal  sphenoidal  array,  the  results  are  not 
complete,  and  since  the  representative  volume  element  of  this  array  is 
not  symmetric,  completion  of  the  stress  strain  relation  in  Makhlouf  and 
Stewart  (1967)  appears  to  be  a  major  task.  The  inexistence  of  the 
primary  source  for  this  array,  Makhlouf  (1963),  made  it  impossible  for 
the  authors  to  clarify  the  above  aspects;  therefore  no  results  will  be 
used  here  for  the  cubical-tetrahedral  and  tetragonal-sphenoidal  array. 


zero,  it  is  low  enough  to  make  this  "vs  =  0  assumption",  needed  for  iso¬ 
tropy,  a  reasonable  one  for  quartz  sands,  at  least  as  a  first  approxima¬ 
tion.  If  vs  =  0  is  assumed,  the  Poisson's  Ratio  of  the  array  is  also 
computed  to  be  v  =  0  for  the  same  three  arrays.  It  must  be  note  that  for 
the  range  0  _<  vs  0.5,  values  up  to  v  -  0.13  are  computed  for  the  same 
arrays  (see  Fig.  32).  Therefore,  the  fact  that  a  value  v  =  0  results  for 
the  array  as  soon  as  vs  =  0  is  assumed  does  not  seem  to  be  so  far  off 
either.  It  is  interesting  to  note  that  measurements  of  Vp  and  Vs  by 
Stokoe  and  his  coworkers  on  an  actual  sand  consolidated  isotropically 
provided  a  similarly  low  value  of  v  =  0.10  (Knox  et  al.  1982,  Kopperman 
et  al.  lq82,  Lee  1985).  In  any  case,  even  with  vs  *  0,  the  cubic 
anisotropy  of  these  expressions  for  v  arrays  is  not  pronounced;  the  error 
resulting  from  computing  the  moduli  between  the  extreme  values  of  vs  does 
not  exceed  3.3%  (Duffy,  1959). 

The  above  three  cubic  arrays,  starting  from  an  isotropic  aQ  state, 
were  loaded  statically  in  triaxial  compression  or  pure  shear  up  to 
failure,  that  is  up  to  gross  sliding,  and  computations  were  performed  and 
are  displayed  here  for  their  stress-strain  curves  and  threshold/failure 
strains  (Figs.  27,  28  and  31).  A  graph  of  obliquity,  022^0  at  failure 
versus  the  intergranular  friction  coefficient  between  spheres,  f,  was 
also  computed  and  is  plotted  in  Fig.  33.  The  curves  obtained  for  the 
arrays  in  this  figure  were  also  compared  to  the  obliquity  obtained 
assuming  the  Mohr-Coulomb  Failure  law: 
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111  =  =  tan2(45  +  4>»  with  tan<f>  =  f 

a0  l-sin$  2 

To  be  able  to  fail  the  simple  cubic  array  in  triaxial  compression,  this 
medium  was  compressed  by  a  force  parallel  to  one  of  the  face  diagonals  of 
the  unit  volume  of  the  array,  that  is  along  a  [110]  direction 
(Deresiewicz ,  1959). 

The  same  cubic  arrays  also  give  excellent  results  when  predicting 
the  influence  of  anisotropic  consolidation  on  shear  wave  velocity;  this 
is  shown  in  Fig.  34  by  a  plot  of  normalized  shear  wave  velocity  vs  stress 
ratio  K  =  022/^0*  In  this  plot,  VS(K)  and  Vs ( 1 )  are  the  values  of  the 
shear  wave  velocity,  Vs,  computed  for  the  anisotropic  case  (K)  and  for 
the  isotropic  loading  condition  (K=l),  respectively,  for  directions  of 
propagation  and  polarization  parallel  or  perpendicular  to  022*  The  same 
plot  includes  data  measured  by  Stokoe  et  al.  (1985)  and  Lee  (1985)  on  dry 
sand  in  the  large  cubic  testing  facility  at  the  University  of  Texas,  with 
excellent  agreement  between  analytical  predictions  and  the  experimental 
data. 

The  shear  modulus  at  very  small  strains,  G^x,  computed  for  these 
same  three  cubic  arrays  under  a  given  isotropic  stress,  oQ,  is  plotted  in 
Fig.  35  as  a  function  of  the  coordination  number  (=  number  of  contacts/ 
sphere).  As  expected,  the  higher  the  coordination  number,  the  stiffer 
the  array,  with  essentially  a  linear  relation  between  the  two  parameters; 
it  is  interesting  that  for  a  given  o0  the  straight  lines  in  Fig.  35  ex¬ 
trapolate  down  to  zero,  suggesting  that  C^ax  is  essentially  proportional 
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to  the  coordination  number  (a  similar  plot  was  obtained  by  Yanagisawa, 
1983).  Therefore,  adding  contacts  to  the  spheres  has  the  same  effect  on 
the  stiffness  of  a  regular  array  as  increasing  the  number  of  springs  in  a 
system  of  equal,  parallel,  elastic  springs.  A  derived  plot  is  the  graph 
between  the  shear  wave  velocity,  Vs  =  (Gmax/p)l/2  versus  void  ratio  and 
isotropic  pressure  o0  (Fig.  36a)  for  the  same  three  arrays.  This  last 
figure  is  especially  interesting,  as  the  trend  observed  in  actual  sands 
is  very  similar  (compare  analytical  curves  with  measured  data  in  sands  in 
Fig.  36b),  except  that  the  absolute  values  of  Vs  in  the  real  soils  are 
much  smaller,  by  a  factor  of  two  or  three.  For  example,  for  e  =  0.47, 
corresponding  to  the  bcc  array,  and  aQ  =  30  psi  =  4,320  psf,  Vs  =  1,800 
fps  is  predicted  by  the  analytical  model  in  Fig.  36a),  while  Vs  =  1,100 
fps  has  been  measured  in  rounded  grained  sands.  Therefore,  Figs.  35  and 
36  strongly  suggest  that  the  dependency  of  Gmax  and  Vs  on  void  ratio 
observed  in  real  soils  is  explained  mainly  by  the  increase  in  the  number 
of  contacts  as  the  void  ratio  decreases. 

Even  though  the  above  results  are  encouraging,  the  regular  arrays 
are  still  very  crude  analytical  models  of  actual  granular  soils,  and 
results  such  as  shown  in  Fig.  36a  are  not  easy  to  interpolate  to  inter¬ 
mediate  void  ratios.  A  significant  improved  model  is  discussed  in  th.‘ 
following  section. 


CHAPTER  6 


AN  ANALYTICAL  MODEL  OF  GRANULAR  SOIL  OF 
ARBITRARY  VOID  RATIO  UNDER  ISOTROPIC  PRESSURE 

Smith  et  al.  (1929)  found  experimentally  that  a  random  arrangement 
of  equal  spheres,  after  enough  shaking  and  tapping  has  been  applied  to 
it,  seems  to  be  composed  of  regular  arrays,  representing  dense  and  loose 
clusters  distributed  within  the  random  grandular  medium.  The  measure¬ 
ments  showed  that  all  spheres  had  between  6  and  12  contacts  per  sphere, 
which  corresponds  exactly  to  the  theoretical  range  for  regular  arrays. 

Additional  experimental  work  by  Bernal  and  Mason  (1960),  Bernal  et 
al.  (1964),  Scott  (1960),  Davis  (1974),  as  well  as  Shahinpoor  and 
Shahrpass  (1982),  and  Finney  (1983),  Figs.  37  and  38),  has  confirmed  that 
2-D  and  3-D  random  assemblages  of  equal  spheres  tend  to  crystalize. 
Consequently,  at  the  present  time,  it  is  generally  accepted  that  an 
assemblage  of  equal  spheres  can  be  modelled  by  a  combination  of  regular 
arrays,  (Finney  1983,  Backman  et  aL.  1983). 

In  this  section,  a  model  of  granular  soil  is  proposed  which  consists 
of  clusters  of  the  three  cubic  arrays  discussed  in  Section  5.5,  with  the 
additional  assumption  that  for  the  spheres  vs  =  0.  In  this  model,  the 
three  cubic  arrays,  having  different  void  ratios  and  inherent  stiffnesses, 
occur  in  proportions  such  as  to  give  the  desired  macroscopic  void  ratio 
of  the  "soil".  In  Section  6.3,  the  relation  between  Graax,  void  ratio  and 
a0  applied  to  this  soil  model  is  calculated  using  the  self  consistent 
method,  and  the  results  are  compared  with  Gmax  measured  in  actual  sands. 
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6.1  The  Self  Consistent  Method 


One  of  the  most  commonly  used  procedure  for  describing  the  behavior 
of  macroscopically  isotropic  composite  elastic  media  is  the  "Self 
Consistent  Scheme". 

This  'Self  Consistent  Scheme"  was  first  devised  by  Hershey  (1954) 
and  Kroner  (1958)  as  a  means  to  model  the  behavior  of  isotropic  and  an¬ 
isotropic  polycrystalline  materials.  Such  materials  are  just  one  phase 
media,  but  because  of  the  random  or  partially  random  orientation  of  the 
crystals,  are  heterogeneous,  their  elastic  properties  vary  with  position 
within  the  medium,  and  discontinuities  in  properties  exist  across  some 
crystal  interfaces. 

In  these  original  applications  of  the  method  to  polycrystalline 
aggregates,  a  single  anisotropic  crystal  was  viewed  as  a  spherical  or 
ellipsoidal  inclusion  within  an  infinite  medium;  this  infinite  medium  had 
the  (still  unknown)  isotropic  elastic  properties  of  the  aggregate.  Then 
the  medium,  with  the  inclusion  in  it,  was  subjected  to  a  uniform  stress 
or  strain  field  applied  at  large  distances  from  the  inclusion.  Next,  the 
orientation  average  of  the  stress  or  strain  inside  the  inclusion  was 
assumed  to  be  equal  to  ("consistent  with")  the  corresponding  applied 
value  of  stress  or  strain.  Thus  the  "self-consistent"  name  of  the 
method.  This  formulation  provided  enough  equations  to  solve  for  the 
isotropic  effective  properties  of  the  medium  (Christensen,  1979). 

Improvement  of  this  self  consistent  scheme  and  its  extension  to 
multiphase  media  are  due  to  Hill  (1965)  and  Budiansky  (1965),  who 
developed  the  method  to  be  used  here.  This  improved  method  represents  an 
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approximate  analysis  for  the  prediction  of  the  overall  (macroscopic) 
elastic  moduli  of  a  multiphase  medium  composed  of  a  coherent  mixture  of 
several  isotropic,  linearly  elastic  materials  or  phases.  The  medium  is 
assumed  to  consist  of  contiguous,  irregular  zones  containing  these 
constituent  materials,  and  the  shapes  of  these  zones  are  assumed  not  to 
deviate  much  from  spherical.  The  spatial  distribution  of  the  phases  is 
assumed  to  be  such  that  the  composite  medium  is  macroscopically  (i.e.  at 
a  scale  much  larger  than  the  dimensions  of  the  zones),  both  homogeneous 
and  isotropic.  Now,  if  an  N-phase  medium  of  total  volume  V  is  defined, 
such  that  the  aggregate  volume  of  all  zones  containing  the  ich  phase  is 
V^,  the  volume  concentration  is  cj  =  V^/V,  and  c^  is  also  equal  to  the 
probability  that  any  arbitrary  point  within  the  medium  ia  located  within 
a  zone  of  the  ictl  material.  It  should  be  noted  that  in  the  limiting  case 
of  very  small  concentrations,  cj,  C2,  •••  cpj-i,  the  first  N-l  phases  will 
tend  to  appear  as  isolated  inclusions  in  a  matrix  consisting  of  the  Nth 
phase. 

In  order  to  obtain  to  effective  overall  (macroscopic)  shear,  G*,  and 
bulk,  K  ,  moduli  of  the  medium,  a  uniform  stress  field  is  applied  at  its 
boundaries.  Then,  the  stress  and  strain  field,  in  each  of  the  phases  is 
evaluated  as  explained  in  the  next  paragraph.  Once  the  fields  are  deter¬ 
mined  for  all  materials,  the  effective  moduli,  G*  and  K*,  can  be  calcu¬ 
lated  by  equating  the  strain  energies  of  the  macroscopic  medium  and  of 
the  phases.  Again,  the  problem  reduces  to  a  number  of  coupled  equations 
tor  K  and  G  ,  which  are  in  terras  of  the  properties  of  the  individual 
materials  and  of  their  volume  concentrations  (Budiansky,  1965).  This 


method  has  been  severely  critized  for  taking  enormous  liberties  with  the 
geometrical  arrangement  of  the  phases  (Christensen,  1979).  To  calculate 
the  elastic  field  in  each  material,  the  geometry  of  the  different  zones 
containing  the  phase  is  successively  rearranged  to  view  the  phase  as  a 
single  inclusion.  However,  the  method  is  relatively  simple  and  in  many 
instances,  when  used  with  caution,  gives  very  good  results.  Furthermore, 
it  has  been  proven  by  Hill  (1965)  that  this  "Self  Consistent  Method" 
yields  results  for  G*  and  K*  which  always  lie  between  the  Voigt  and  Reuss 
bounds,  that  is,  the  spatial  average  of  the  moduli  of  the  phases  (Voigt 
bound,  springs-in-parallel)  and  of  the  reciprocal  of  the  moduli,  or 
compliances  of  the  phases  (Reuss  bound,  springs-in-series) . 

The  evaluation  of  the  stress  and  strain  fields  in  each  of  the  phases 
is  performed  for  the  isotropic  case  by  the  solution  of  the  problem  of  the 
elastic  field  inside  an  ellipsoidal  elastic  inclusion  (Eshelby,  1957). 

It  was  shown  by  Eshelby  that  the  elastic  field  inside  an  ellipsoidal 
isotropic  elastic  inclusion  embedded  in  an  isotropic  elastic  medium  is 
uniform;  this  is  an  extremely  important  conclusion  as  it  eliminates  the 
need  for  averaging  the  fields  within  the  inclusion  phase  and  simplifies 
enormously  the  formulation.  Later,  it  was  shown  that  the  stress  and 
strain  fields  inside  an  orthotropic  inclusion  embedded  in  an  orthotropic 
medium  are  also  uniform,  as  long  as  the  cross  section  of  the  inclusion  is 
quadratic  (Kinoshita  and  Mura,  1971). 

The  averaging  of  the  shear  moduli  of  all  phases  by  a  strain  energy 
balance  between  the  medium  and  the  inclusions  yields: 
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N— 1 


1.1.  V  n 

m  +  Y  Ci[l  -  —  )(-o) 
gN  i=  1  %  T 


(42) 


i  .  N_1  k  .  £v, 

p-f +  i  <=i(‘ 

K  i=l  KN  a0 


(43) 


where  G*,  K*  are  the  desired  macroscopic  moduli;  Kj,  (i  =  1,2, ...N) 

are  the  moduli  of  the  ic^  phase,  Cj_  =  V^/V  is  the  volume  concentration, 

and  Y£,  ev^  are  the  values  of  the  average  shear  strain  and  volumetric 

strain  respectively,  inside  the  phase.  The  parameters  t°  and  o°  are  the 

o 

shear  stress  and  isotropic  pressure  applied  at  the  boundary  of  the 
medium. 

The  Elshelby  (1957)  solution  gives 


Yi 


G*+0*(Gi-G*) 


(44) 


1  K*+a*(Ki-K*) 


(45) 


ft  If 

where  ot  ,8  are  components  of  the  Eshelby  S-tensor;  for  the  case  of 
spherical  inclusions  they  are: 


l+V 


3( 1-v*) 


(46) 


ft*  a  2  (4-5v  ) 
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(47) 


where  v  is  the  macroscopic  Poisson's  ratio  of  the  medium: 


By  "smearing  out",  that  is,  by  replacing  the  matrix  surrounding  each 
inclusion  (phase)  by  the  desired  resultant  macroscopic  medium,  equations 
(42)  and  (43)  reduce  to: 


N  c 

l  - ~ 

i="l  *rGi 


y  ci 

i -1  n 

1  +a*  ( 0 


which  are  symmetrical  for  the  various  phases.  Therefore,  Budiansxy  (1965) 
has  suggested  to  use  Equations  49  and  50  for  arbitrary  concentrations  of 
the  constituents  of  the  composite  medium  as  described  previously.  Fur¬ 
thermore,  Budiansky  simplified  equations  44  and  45  to: 


r - i - 

G*  1 

iV(^“  0 


1+0  K* 


A  comparison  of  Eq.  51  with  results  obtained  through  statistical 
finite  element  methods,  suggests  that  the  above  equations  do  indeed  model 
the  continuum  described  previously,  including  the  assumption  that  the 
stress  and  strain  fields  of  the  phases  are  approximately  independent  of 


n 

$ 


35* 


i 


m 

sm 


location  (Fig.  39,  see  also  Petrakis,  1983). 

Finally,  Equations  46-50  have  to  be  solved  simultaneously  to  yield 
the  desired  values  of  the  macroscopic  elastic  moduli,  G*,  K*.  These 
resultant  macroscopic  moduli  are  estimates  of  the  overall  elastic 
constants  of  the  multiphase  medium,  and,  as  mentioned  before,  they 
invariably  lie  between  the  Reuss  and  Voigt  bounds.  Other  solutions  may 
provide  narrower  bounds  for  the  actual  solution  (Hashin  and  Shtrikman, 
1963);  however,  the  Self  Consistent  solution,  in  certain  cases,  also 
falls  between  these  narrower  bounds,  thus  showing  its  capability  for 
providing  accurate  results  (Hill,  1965). 

6 . 2  The  Model 

The  Self  Consistent  Method  is  applied  here  to  evaluate  the  elastic 
constants  of  a  random  assemblage  of  equal,  rough  elastic  spheres  that  has 
been  consolidated  isotropically  and  has  a  prescribed  mean  void  ratio  e. 
The  spheres  are  assigned  the  elastic  properties  of  quartz,  and  the  assem¬ 
blage  is  assumed  to  be  composed  of  random  zones,  with  each  zone  consist- 
uted  by  a  large  number  of  spheres  arranged  in  either  of  three  regular 
cubic  arrays. 

Recently,  Shahinpoor  (1981)  modelled  a  random  2-D  array  of  equal 
steel  spheres  as  a  combination  of  Voronoi  cells,  derived  an  expression 
for  the  probability  density  function  of  the  void  ratio,  p(e),  and  checked 
experimentally  this  analytical  p(e)  by  means  of  an  optical  scanning 
technique,  Fig.  37  (see  also  Shahinpoor  and  Shahrpass,  1982).  The 
expression  for  p(e)  is: 


,’V,4  -  r'j.l't.f  .<  .tl  >  *“•  |k*.  f.S  |#»  |»*  h*  l.t  l.|  d.  I  ft  I'.  |l.  IH  *ft»  1.1  , 


ft 


p(e) _ X  exp(-Xe) _ 

exp(-Xem£n)-exp(-Xemax) 


.  -  1  emin  expC-Xeminj-emax  exp(~Xemax) 

w  ere  e  -  A  +  exp( -Xeraln)-exp( -Xeraax) 


ft 

$ 

$ 

v! 


X  is  obtained  from  the  mean  void  ratio,  e,  of  the  distribution  p(e).  As 
mentioned  before,  it  is  reasonable  to  model  a  uniform,  rounded-grained 


sand  as  a  random  combination  of  zones  corresponding  to  regular  cubic 
arrays,  and  this  was  the  approach  taken  in  this  work.  The  sand  medium  is 
assumed  to  be  composed  of  regular  arrays  in  the  fashion  of  Figs.  37  and 


38,  where  each  randomly  oriented  Voronoi  polyhedron  is  one  of  these 
zones,  and  contains  a  regular  array  with  many  spheres.  A  cross  section 
of  this  3-D  medium  could  be  visualized  approximately  by  the  actual 
photograph  of  the  2-D  medium  in  Fig.  37b;  in  this,  the  black  spots  are 
spheres  and  the  white  are  voids,  and  zones  of  regular  packings  can  be 
clearly  observed.  The  macroscopic  moduli  of  the  whole  medium  will  be 


determined  from  the  properties  of  these  zones  through  the  self  consistent 
s  cheme . 

As  a  first  step,  the  probability  density  function  of  the  void  ratio, 
p(e),  Eq.  53,  was  transformed  to  the  probability  density  function  of  the 
porosity,  p(n),  with  the  basic  equation  (Benjamin  and  Cornell,  1970): 

P(n)  =  jff |  P(e ) 

as  the  mean  of  the  porosity  distribution,  n,  coincides  with  the  macro¬ 
scopic  (measured)  porosity  of  the  "soil",  unlike  the  mean  void  ratio  e, 
which  is  not  identical  to  the  macroscopic  void  ratio  (Petrakis,  1983). 


Then,  the  probability  density  function  of  the  porosity,  p(n),  was 
discretized  into  three  segments  of  influence,  corresponding  respectively 
to  the  porosities  of  the  three  regular  cubic  arrays  (see  Fig.  40):  sc 
(n  =  0.48),  bcc  (n  =  0.32)  and  fee  (n  =  0.26).  The  values  of  nmj[n  and 
nmax  used  for  all  calculations  were  those  of  the  sc  and  fee  arrays. 

For  example,  for  a  prescribed  macroscopic  porosity  n  =  0.35  corres¬ 
ponding  to  a  mean  void  ratio  e  =  0.54,  the  calculation  illustrated  in 
Fig.  40  allowed  determining  the  following  three  volume  concentrations, 
ci : 


Array 

^i_ 

sc 

0.48 

0.1934 

bcc 

0.32 

0.5921 

fee 

0.26 

0.2143 

The  medium  with  these  three  phases  was  then  subjected  to  an 

isotropic  boundary  confining  pressure,  a°,  and  subsequently  subjected  to 

o 

small  boundary  stress  increments  da^j,  from  which  the  corresponding 
elastic,  very  small  strain  increments  at  the  boundaries,  and  the  mean 
cubic  medium  moduli  K*  and  G*  were  evaluated. 

If  we  now  assume  that  the  phases  are  quadratic  (elliptical  or 
circular)  in  cross  section  we  can  apply  the  Self  Consistent  method.  The 
assumption  that  the  "zones"  are  quadratic  in  cross-section  is  important, 
since  if  the  resulting  stress  and  strain  fields  are  uniform  (see  Section 
6.1)  within  each  "zone”,  the  "zone"  can  be  replaced  by  the  representative 
cube  of  each  cubic  array  (Figs.  24-26)  and  the  corresponding  constitutive 
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laws  are  given  by  Eqs.  18,  24  and  31.  Since  in  turn  these  relations 
depend  upon  the  pressure  acting  on  each  inclusion  or  phase,  the  value  of 
the  stress  field  at  the  boundary  of  each  of  these  phases  and  inside  it 
can  be  readily  obtained  from  Eshelby's  (1957),  and  Budiansky's  (1965) 
results: 


-i 

<*o 


+  a*(^  '  1) 


(55) 


where  the  value  of  a*  depends  upon  the  shape  of  the  zone,  (Eshelby, 

1957),  which  here  has  been  assumed  to  be  spherical  for  simplicity.  Note 
_i 

that  this  value  of  oQ  is  independent  of  the  location  of  the  zone,  and  is 
thus  the  same  for  all  zones  containing  the  same  regular  cubic  array  or 
phase. 

Equation  55  is  then  replaced  into  Eqs.  18,  24  and  31,  and  the 
problem  finally  reduces  to  the  solution  of  the  following  equations  for 
the  three  phases: 


_i 

Oo 


o  Kt 

°o  [- 


1 


1  + 


Kt 


J  i  =  1,2,3 
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=  1.0 


(58) 


-  1) 


1+v 


3( 1-v  ) 


(59) 


8  =  TT 


2  (4-5v*) 

15  (1-v*) 


(60) 


3K*-2G* 


6K*+2G* 


(61) 


where  Kj  =  ^(oq),  G^  =  G£(a0)  for  the  three  phases  are  obtained  with  Eqs. 
10,  16  and  23  for  i  =  1,  2,  3,  corresponding,  respectively,  to  the  simple 
cubic,  body  centered  cubic  and  face  centered  cubic  regular  arrays. 


6.3  Application  to  Quartz  Sand 

The  proposed  model  was  evaluated  using  as  input  the  elastic  param¬ 
eters  of  quartz  for  the  individual  spheres,  which  are  Es  =  11,0  x  10^  psi 
and  v3  =  0.15  (White,  1965,  Ko  &  Scott,  1967,  Lambe  and  Whitman,  1969, 
see  Table  3),  and  for  a  wide  range  of  isotropic  confining  pressures.  The 
values  of  the  computed  shear  modulus  G*  are  plotted  on  solid  lines  in 


Figs.  41,  42  and  43  versus  confining  pressure  o0  =  a0,  for  e  =  0.46,  e  = 

o 


0.54,  and  e  =  0.58,  respectively.  The  values  of  the  bulk  modulus,  K  , 
were  also  computed,  and  Fig.  44  contains  a  plot  of  confining  pressure 
versus  volumetric  strain  predicted  by  the  model  for  e  =  0.54,  where  the 


volumetric  strain  was  derived  from  this  computed  bulk  modulus,  ev  =  o0/K 


As  mentioned  before  in  connection  with  Fig.  36,  regular  arrays  are 
much  stiffer,  up  to  3.5  times  stiffer,  than  actual  uniform,  rounded 
sands,  and  this  constitutes  a  serious  defficiency  of  the  model.  However, 
Drnevich  and  Richart  (1970)  have  suceeded  in  increasing  significantly  the 
shear  stiffness  of  dry,  rounded,  uniform  Ottawa  quartz  sand  by  applying 
millions  of  cycles  of  a  shear  strain  slightly  greater  than  the  threshold 
strain  in  a  resonant  column  device. 

Figures  41,  42  and  43  also  include  as  data  points  the  resonant 
column  experimental  results  obtained  by  Drnevich  (1967),  and  Drnevich  and 
Richart  (1970),  for  the  same  macroscopic  void  ratios,  e  =  0.46,  e  =  0.54 
and  e=  0.58,  used  to  calculate  the  solid  lines.  The  lower  dotted  line  in 
each  figure  corresponds  to  virgin  or  uncycled  sand  as  predicted  by  the 
Hardin  and  Black  (1966)  correlation,  which  is  essentially  identical  to 
the  Hardin  and  Richart  correlation  depicted  in  Fig.  36(b).  In  Figs.  41- 
43,  two  trends  may  be  clearly  observed:  a)  as  the  number  of  cycles,  N, 
increases,  the  test  results  steadily  approach  the  model  values  until,  at 
N  =  1  to  2  x  107  cycles  the  agreement  becomes  excellent;  and  b)  the  slope 
of  the  line  of  shear  modulus  vs  confining  pressure  decreases  from  about 
1/2  in  the  uncycled  state  to  the  theoretical  1/3  after,  approximately,  1 
x  10^  cycles.  The  reason  for  some,  of  the  points  showing  more  scatter  is 
probably  because  those  points  were  cycled  more  than  others  with  the  same 
void  ratio,  or  because  their  void  ratios  were  slightly  differed. 

In  the  above  tests,  the  sand  specimens  were  cycled  at  strains  which, 
although  larger  than  the  threshold  strain,  were  small  enough  so  that  no 
significant  densif ication  occurs,  and  indeed  the  change  in  measured 
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(macroscopic)  void  ratios  between  virgin  and  cycled  specimens  was  very 
little  or  negligible;  thus,  densif ication  is  certainly  not  the  explana¬ 
tion  for  the  observed  threefold  increase  in  the  sand  stiffness.  Drnevich 
and  Richart  (1970)  speculated  in  their  paper  that  the  above  behavior 
could  be  due  to  wearing  of  the  contacts,  increase  of  the  contact  areas  or 
formation  of  additional  contacts.  The  authors  think  that  this  third 
reason  explains  the  phenomenom  completely,  as  illustrated  by  the  compar¬ 
ison  with  the  model  in  Figs.  41-43,  and  by  the  relation  between  stiffness 
and  number  of  contacts  in  Fig.  35.  It  is  known  that  in  a  random  array  of 
spheres  there  may  be  less  contacts  than  in  regular  packings  (see  Smith  et 
al.,  1929)  and,  furthermore,  it  is  possible  to  have  contacts  which  do  not 
transmit  any  load  (dead  contacts).  By  continuous  cycling  such  as 
performed  by  Drnevich  and  Richart  in  their  tests  these  contacts  were  made 
load-transmitting  and  new  ones  were  formed  until  all  or  most  possible 
contacts  were  created  and  the  stiffness  of  the  sand  coincided  with  that 
predicted  by  the  model.  It  is  interesting  that  the  creation  of  contacts 
can  also  be  achieved  by  high  isotropic  confining  pressures  (Duffy  and 
Mindlin,  1958,  Deresiewicz  1958,  see  also  Fig.  13). 

The  analytical  model  proposed  herein  describes  exactly  this:  since 
it  assumes  that  the  sand  is  a  random  assemblage  of  regular  arrays  of 
spherical  grains,  it  implies  that  the  number  of  contacts  is  the  maximum 
possible.  Furthermore,  as  shown  by  Duffy  and  Mindlin  (1957)  and 
Deresiewicz  (1958),  ar.d  as  illustrated  by  Fig.  13,  when  the  number  of 
contacts  increases  the  pressure  dependency  of  the  moduli  tends  co  change 
from  1/2  to  1/3.  All  of  this  is  interpreted  by  the  model,  which  is  shown 


4,« 


here  to  represent  the  limiting  state,  in  terms  of  number  of  contacts  and 
of  stiffness,  that  a  sand  can  reach. 

o  _ 

Figure  44  shows  the  confining  pressure  o0  vs.  volumetric  strain,  ev, 
measured  by  Drnevich  and  Richart  on  the  same  Ottawa  sand  specimens  dis¬ 
cussed  above,  for  a  virgin  specimen  and  for  a  specimen  after  1  x  10^ 
cycles  of  shear  strain.  Unfortunately,  no  data  is  available  for  sand 
specimens  cycled  with  more  than  1  x  10^  cycles.  In  the  figure,  lines 
have  been  passed  which  approximately  represent  these  experimental  data. 
The  line  predicted  by  the  model  is  also  plotted,  and  again  it  is  clear 
that  the  cycling  increased  the  bulk  modulus  of  the  soil,  thus  making  it 
approach  the  predicted  model  curve,  as  the  number  of  contacts  increases 
toward  the  theoretical,  maximum  value.  It  would  be  interesting  to 
compare  the  difference  between  the  analytical  results  and  the  experimen¬ 
tal  values  for  both  moduli,  K*  and  G*  at  a  given  cycling  state.  The 
experimental  values  for  the  case  of  hydrostatic  compression  (Fig.  44)  are 
closer  to  the  model  predicted  curve  than  the  measurements  with  shear 
(Figs.  41-43),  for  a  comparable  number  of  cycles.  This  is  explained  by 
the  fact  that,  once  the  spheres  have  approached  each  other  due  to  the 

cycles  of  shearing,  an  increase  in  o°  can  complete  the  formation  of  many 

o 

contacts,  thus  further  increasing  the  stiffness  of  the  soil. 
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CHAPTER  7 

A  TWO-DIMENSIONAL  NUMERICAL  MICROMECHANICAL  MODEL  FOR  THE 
STRESS-STRAIN  BEHAVIOR  OF  A  GRANULAR  SOIL  AT  SMALL  STRAINS 


In  the  preceeding  chapter  a  closed  form  solution  was  obtained  for 
the  elastic  constants  of  a  random  configuration  of  equal,  rough,  elastic 
spheres  having  an  arbitrary  macroscopic  void  ratio  and  subjected  to 
isotropic  loading.  This  was  done  through  the  use  of  the  Self-Consistent 
Method,  using  a  minimum  number  of  assumptions,  and  it  contributed  to  a 
greatly  improved  understanding  of  the  small  strain  behavior  of  sands 
under  isotropic  conditions.  The  fact  that  this  analytical  solution  was 
obtained  with  relatively  small  effort  should  be  attributed  to  the  high 
level  of  symmetry  in  such  a  system  under  isotropic  pressure. 

Logically,  the  next  step  should  have  been  to  investigate  the 
possibility  of  obtaining  equivalent  rigorous  solutions  for  the  case  of  a 
similar  random  configuration  of  equal,  elastic,  rough  spheres,  but  now 
subjected  to  anisotropic  loading.  Unfortunately,  this  proved  impossible, 
not  only  because  of  the  lower  level  of  symmetry  in  this  case,  which  made 
the  combination  of  three  or  more  phases  for  the  medium  impossible,  but 
also  because  anisotropic  loading  induces  early  yielding  (sliding)  between 
some  particles  in  the  individual  phases.  Therefore,  an  elastic  solution 
was  not  possible.  An  elastic,  perfectly-plastic  rigorous  solution  would 
have  been  -in  principle  feasible,  but,  again,  the  pressure  dependence  of 
the  properties  of  the  phases  made  the  problem  intractable. 
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As  mentioned  in  Section  4.2,  for  years  researchers  have  been 
modeling  the  elastoplastic  behavior  of  polycrystalline  aggregates  using  a 
variety  of  analytical,  semianalytical,  and  numerical  techniques.  As 
discussed  previously,  exhibit  an  analogous  behavior  to  that  of  the  sands. 

Therefore,  it  seems  logical  to  attempt  formulating  the  problem  of  a 
random  packing  of  equal  spheres  subjected  to  an  arbitrary  stress  path,  in 

i 

a  manner  similar  to  that  used  for  crystals.  This  can  be  achieved  by 

devising  an  aggregate  of  regular  arrays  such  as  that  of  Chapter  6, 

describing  the  force-deformation  at  the  intergranular  contacts  through 

the  nonlinear  numerical  techniques  described  in  Section  4.1  (subroutine  , 

i 

CONTACT),  and  resorting  to  a  nonlinear  finite  element  simulation. 

7.1  Aggregate  Description 

In  order  to  model  the  behavior  of  granular  soil  at  small  strains,  a  J 

numerical  simulation  was  developed  which  calculates  the  response  of  a 

random  aggregate  of  equal,  elastic,  rough  spheres  under  an  arbitrary 

o 

boundary  stress  state,  o  .  For  this,  a  finite  element  analysis  was  per- 

ij 

formed  in  which  the  element  corresponded  to  a  simple  cubic  array.  Each 
element  contains  an  undetermined  number  of  spheres,  assumed  to  be  subjec¬ 
ted  to  a  uniform  stress  field,  so  that  four  spheres  represent  the  above 

I 

element  (Duffy  and  Mindlin,  1957);  the  element  configuration  appears  in  j 

Fig.  45.  Several  specific  media  were  used  in  these  simulations,  each  ] 

medium  consisting  of  a  number  of  simple  cubic  arrays  in  two  dimensions  (a  , 

mono-layer  of  equal  spheres  assembled  in  simple  cubic  patterns),  oriented  J 

in  such  a  way  so  as  to  resemble  a  statistically  isotropic  random  aggregate. 
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The  effect  of  the  shape  of  the  "grain  boundaries",  in  this  case  the  gap 
between  elements  which  are  oriented  differently,  and  of  any  localized 
orientation  concentrations  on  the  macroscopic  behavior  was  not  studied 
here,  since  a  statistically  uniform  medium  was  sought  to  interpret  the 
macroscopic  behavior  of  a  granular  soil.  Each  of  the  above  elements  is  a 
simple  cubic  array,  randomly  oriented,  and  the  number  of  elements  used  in 
each  of  these  media  varies  from  16  (4X4)  to  64  (8X8)  as  illustrated  in 
Fig.  46  and  Figs.  D1  through  D2  in  Appendix  D.  Therefore,  the  medium  has 
the  same  properties  of  the  simple  cubic  array  in  terms  of  void  ratio  and 
coordination  number. 

The  simple  cubic  array  was  selected  as  the  basis  for  all  simulations 
due  to  its  simplicity.  It  is  extremely  difficult,  for  example,  to  imple¬ 
ment  numerically  a  constitutive  law  for  the  body-centered  cubic  array 
since  the  contact  forces  and  displacements  are  coupled.  For  another 
alternative,  the  face-centered  cubic  array,  the  task  is  impossible  as  the 
array  is  statically  indeterminate.  However,  it  was  felt  that  the 
inability  of  varying  the  void  ratio  or  the  coordination  number  (number  of 
contacts  per  particle)  did  not  unduly  constrain  the  simulations,  which 
were  expected  to  provide  useful  Insights  into  the  behavior  of  sand  at 
small  strains. 


7.2  The  Element 

The  constitutive  relation  of  the  elements  described  above,  is  an 
incrementally  linear  law  which  is  based  on  the  differential  stress-strain 
expressions  for  the  Simple  Cubic  array  presented  in  Chapter  5.  These 
incremental  relations  at  the  particle  contacts  are  integrated  through  the 


nonlinear  program  CONTACT,  previously  described  in  some  detail  in  Section 
4.1.  Thus,  at  every  increment,  the  displacements  in  the  simple  cubic 
array  are  computed  separately  for  each  contact,  and  then  the  results  are 
combined  following  the  formulation  presented  in  Chapter  5. 

Since  finite  element  programs  are  implemented  through  the  Stiffness 
Method,  in  which  the  displacements  are  the  independent  variables,  and 
thus  the  code  first  computes  the  displacements  at  the  nodes  from  the 
specified  boundary  conditions  and  then  calculates  the  stresses  at  the 
Gauss  points  from  these  nodal  displacements,  the  constitutive  relation  of 
the  incremental  model  CONTACT  (Eqn.  12)  had  to  be  inverted  in  order  for 
the  displacements  to  become  the  independent  variables.  The  inverted 
constitutive  relation  is  as  follows: 


+  JCji  f  HN 

dF  -  da  k  +  Hp  dDt  t  +  fdN  n  +  H(dDN  -  -g22-)  n 


if  dT/dN  >  f 


+  2Ga 

dF  -  da  k  +  Hp  dDc  t  +  H0  dD^  n 


if  dT/dN  <  f. 

In  the  general  case  in  which  the  simple  cubic  array  element  is  com¬ 
pressed  biaxially  under  aj[  and  a 22  (see  Fig.  45),  the  normal  strains  are 
uncoupled  and  depend  only  on  the  normal  relative  displacements  between 
the  centers  of  adjacent  spheres,  which  in  turn  depend  only  on  the  normal 
stresses.  As  described  in  Chapter  4,  this  occurs  because  the  behavior  at 
a  contact  subjected  to  a  normal  force  is  nonlinear  elastic,  and  because 
the  normal  compliance  at  the  contact  is  only  functions  of  the  normal 
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force.  Unfortunately ,  this  Is  not  true  for  the  shear  strains,  which  are 
coupled  with  the  normal  stresses.  The  expression  for  the  shear  strain, 
Y12»  described  in  Section  5.1,  is  yj2  *  C 6 1 2  +  <$2l)/2R,  where  612  and  62 1 
are  the  relative  tangential  displacements  between  the  centers  of  adjacent 
spheres  of  radius  R.  Therefore,  the  tangential  displacements,  612  an<* 
$21»  depend  on  the  values  of  both  the  tangential  and  the  normal  forces  at 
the  contact. 

As  a  first  step,  the  macroscopic  normal  stresses  corresponding  to 
each  contact  are  computed  directly,  since  the  behavior  of  the  array  under 
compression  is  nonlinear  elastic.  Once  these  normal  stresses  and  forces 


are  determined,  the  tangential  forces  and,  consequently,  the  shear 
stresses  can  also  be  computed.  The  shear  strain,  Y12,  supplied  by  the 
finite  element  program  at  each  Gauss  point,  is  obtained  from  the 
summation  of  the  two  tangential  displacement  components  at  both  contacts, 
6 12  ar*d  ^2l«  Since  the  shear  stress  components  on  the  faces  of  the 
element,  012*  must  be  equal,  the  problem  of  finding  the  value  of  012 


reduces  to  computing  the  common  tangential  force  acting  on  two  different 
contacts,  subjected  to  two  different  normal  forces,  and  with  the 
summation  of  the  resulting  tangential  displacements  being  known  (<Sj2  + 


621  3  2Ryi2)*  For  this  purpose,  an  iterative  numerical  technique  based 
on  the  bisection  method  was  used,  and  this  made  possible  the  evaluation 
of  the  tangential  force  and,  subsequently,  of  the  shear  stress. 

The  stress-strain  behavior  of  an  individual  element  (Fig.  45) 
assuming  quartz  spheres,  is  shown  in  Fig.  47c  for  a  biaxial  state  of 
stress,  on  and  022*  followed  by  pure  shear.  Figures  47a  and  47b  portray 
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the  force-deformation  behavior  of  the  "weak"  and  "strong"  contacts  for 
the  same  case,  respectively,  where  the  "weak"  contact  is  defined  as  that 
for  which  the  normal  stress  is  smaller  (on  in  this  case).  The  proper¬ 
ties  of  quartz  for  the  spheres  utilized  in  this  calculation  (E  =  295,182 
Kg*/cm2,  f  =  0.5  and  v  =  0.15,  White,  1964)  are  consistently  used 


throughout  the  remainder  of  this  work. 

It  should  be  noted  that  the  element  can  not  take  any  tensile  normal 
stress  on  the  slip  planes,  since  this  would  imply  that  a  contact  has 
ceased  to  exist,  and  thus  that  the  particles  want  to  rearrange  themselves 
and  form  a  new  packing*.  Although  this  rearrangement  obviously  happens 
in  actual  sand  aggregates,  its  simulation  was  too  complex  to  implement  in 
the  present  finite  element  code,  and  the  decision  was  made  to  make  the 
contact  normal  forces  positive-definite,  that  is  they  can  be  positive  or 
zero,  but  never  negative. 


A  subroutine  implementing  the  behavior  of  this  element  was  coded 
into  the  nonlinear  finite  element  program  ABAQUS  (1982)  as  a  user  defined 
material  subroutine  (UMAT).  The  listing  of  the  subroutine  as  well  as  a 
flowchart  appear  in  Appendix  C.  Finally,  an  incrementally  linear 
analysis  was  done  using  plane  strain  eight-noded  elements  with  reduced 
integration. 


*For  the  simple  cubic  array  element  considered  herein,  simple  analytical 
considerations  show  that  if  both  oi[  and  <722  are  positive  definite,  then 
no  tensile  stress  can  exist  in  any  direction  if  the  friction  coefficient 
f  is  less  than  1,  as  is  the  case  here. 
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CHAPTER  8 

MONOTONIC  LOADING  SIMULATIONS  OF  THE  AGGREGATE 

As  a  first  step,  a  number  of  monotonic  loading  cases  were  simulated 
in  an  attempt  to  define  the  macroscopic  behavior  of  the  aggregate.  The 
medium  was  subjected  to  several  loading  paths  ranging  from  hydrostatic 
compression,  to  hydrostatic  compression  followed  by  some  combination  of 
shear  and  compression. 

In  most  of  these  simulations,  an  attempt  was  made  to  keep  the  mean 
stress  constant,  so  as  to  have  the  simulations  contained  on  one  x-plane, 
and  the  microscopic  behavior  was  closely  monitored  in  order  to  verify 
that  the  necessary  assumptions  of  the  constitutive  law  were  satisfied. 
Finally,  one  simulation  was  also  made  with  variable  mean  stress  in  order 
to  investigate  the  behavior  of  the  aggregate  under  biaxial  compression. 

8 . 1  Isotropic  Compression 

In  one  of  the  first  simulations,  medium  1  with  16  elements  was  sub¬ 
jected  to  a  monotonically  increasing  hydrostatic  pressure,  first  to 

°  ?  °  9 

oQ  «  3  Kg*/cmz  and  then  up  to  o0  ■  5  Kg*/cnr,  and  it  was  observed  that 
both  the  macroscopic  and  microscopic  (element)  response  exhibit  a  locking 
nonlinear  elastic  behavior,  similar  to  that  observed  in  sand  (Fig.  48). 

8.2  Isotropic  Compression  Followed  by  Pure  Shear 

As  a  second  step,  the  isotropy  of  the  media  used  was  verified  by 

o  9 

loading  them  isotropically  up  to  o0  ■  1  Kg*/cmz,  followed  by  pure  shear 
applied  incrementally.  This  was  accomplished  by  imposing  on  each  medium 
a  predetermined  direction  of  the  major  principal  stress.  The  values  used 
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for  the  angle,  a,  between  the  major  principal  stress  and  the  vertical 

direction  of  the  medium,  (Fig.  49a),  were  0°  (compression  in  the  vertical 

direction);  22.5°;  45°  (pure  shear  in  the  vertical  and  horizontal 

planes);  67.5°;  and  90°  (extension  in  the  vertical  direction).  The 

results  of  these  simulations  are  shown  in  Fig.  50  as  plots  of  the  applied 
o  o 

deviator  stress,  04-02,  versus  the  resultant  shear  strain,  y  =  ej-e2» 

Fig.  50  includes  results  of  media  with  16  and  64  elements.  In  total, 
four  media  were  used:  three  with  16  and  one  with  64  elements,  including 
the  two  media  in  Fig.  46.  Figs.  51-53  include  the  stress-strain  behavior 
of  two  individual  elements  for  the  compression  of  medium  2  under  constant 
mean  stress  and  for  a  -  0°.  The  stress-strain  behavior  of  all  16  ele¬ 
ments  for  the  same  medium  and  same  loading  case  appears  in  Appendix  E. 

It  can  be  seen  in  Fig.  50  that  the  aggregate  is  indeed  isotropic 
under  isotropic  pressure,  as  expected.  Since  Fig.  50  shows  that  the 
difference  between  the  stress-strain  behavior  of  the  16-element  and  64- 
element  media  is  not  significant,  it  was  decided  that  for  subsequent 
parametric  studies,  as  well  as  for  monitoring  the  stress-strain  behavior 
of  each  element,  any  of  the  less  costly  16-e lenient  media  could  be  used  as 
representative  of  the  aggregate. 

6.3  Isotropic  Compression  Followed  by  Biaxial  Compression 
A  16-element  medium  (medium  1  of  Fig.  46)  was  compressed  with  the 
mean  stress  variable,  in  a  manner  similar  to  that  of  a  soil  specimen  in  a 
triaxial  device  (see  the  45°  stress  path  in  Fig.  54).  That  is,  the 

o  9 

medium  was  first  compressed  isotropically  under  aQ  =  1.0  Kg*/cmz,  and 


then  oi  was  increased  in  a  vertical  direction  while  02  remained  constant, 
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°2  “  °o  a  1*0  Kg*/cmz.  The  stress-strain  curve  calculated  for  this 
medium  is  shown  in  Fig.  55,  and  the  behavior  for  two  individual  elements 
is  presented  in  Figs.  56  and  57.  The  a~e  behavior  of  all  16  elements 
for  the  same  medium  and  same  loading  case  appears  in  Appendix  F.  The 
nonlinearity  of  the  medium  here  is  even  less  pronounced  than  in  the 
simulations  with  constant  mean  stress  (Fig.  50),  since  in  Fig.  55  there 
is  no  lateral  unloading.  The  Poisson's  ratio  of  the  medium  was  also 
computed  from  this  numerical  simulation,  and  it  was  found  to  vary  from 
0.02  to  0.05  (Fig.  58),  thus  agreeing  well  with  values  measured  experi¬ 
mentally  on  dry  sand  at  the  University  of  Texas  (Lee,  S.,  1985).  The 
obliquity  at  failure  was  the  same  as  in  the  simulations  with  the  constant 
mean  stress,  as  illustrated  in  Fig.  54. 

8.4  Yielding  and  Failure  Considerations 
In  the  numerical  simulations  of  pure  shear  summarized  in  Fig.  50,  it 
was  observed  that  the  yielding/failure  process  of  the  medium  occurs  in 
two  successive  stages.  In  the  first  stage,  a  growing  number  of  "soft" 
elements,  oriented  more  or  less  parallel  to  the  directions  of  the  applied 
shear  stress,  slide  and  this  sliding  accounts  for  the  increasing  nonlin¬ 
earity  of  the  curves  in  Fig.  50,  as  the  shear  strain,  y  *  ei~e2»  increases 
from  0  to  values  around  0.1  x  10-3.  At  these  larger  values  of  shear 
strain,  typically  around  20%  of  the  elements  have  already  slid  (failed). 

In  the  second  stage,  occurring  at  about  y  =  ei-e2  =  0.16  x  lO--^,  one  or 
several  of  the  "stiff"  elements,  oriented  more  or  less  perpendicular  to 
the  direction  of  the  shear  stress,  and  which  had  not  yet  slid,  tend  to 
separate  as  the  normal  force  at  the  contact  becomes  zero,  the  correspond- 


ing  ratio  shear/normal  force  at  the  contact  reaches  f,  and  the  element 
slides.  This,  of  course,  is  related  to  the  fact  that  the  normal  contact 
force  is  allowed  to  be  zero  but  not  negative.  Once  some  of  the  "stiff" 
elements  fail  due  to  this  tendency  to  separate,  a  growing  number  of  both 
"soft"  and  "stiff"  elements  slide  in  the  next  increment(s)  by  a  combina¬ 
tion  of  shear  stress  increase  and  separation  tendency,  thus  precipitating 
the  failure  of  the  medium,  occurring  at  if  =  0.448  Kg*/m2.  This  is 
close,  but  slightly  less  than  the  yield  stress  of  the  simple  cubic  array 
subjected  to  pure  shear  in  the  directions  of  the  array:  t  =  (0.5)/ (1)  = 
0.5  Kg*/cm2.  This  value  of  the  failure  stress  is  much  less  than  the 
value  of  1.536oy  calculated  for  the  polycrytalline  aggregates,  as 
discussed  in  Section  4.2.  This  difference  in  values,  between  the  simple 
cubic  and  the  fee  polycrystalline  aggregate,  occurs,  among  other  reasons, 
because  the  simple  cubic  element  is  not  allowed  to  take  any  tensile 
stresses  during  the  stress  redistribution  which  takes  place  when  an 
increasing  number  of  the  elements  slide. 

To  further  verify  this  sequence  of  events  just  described,  a  special 
64-element  medium  was  defined  (Fig.  D3  in  Appendix  D)  and  subjected  to 
shear  at  constant  mean  stress.  This  new  medium  was  constructed  in  such  a 
way  as  to  provide  an  insight  into  this  complex  "frilure"  phenomenon,  by 
having  only  two  extreme  element  orientations:  the  four  elements  at  the 
center  are  oriented  at  8  =  0°,  see  Fig.  49b  ("stiff"  elements),  while  the 
other  sixty,  surrounding  those  four,  are  oriented  at  45°  ("soft"  ele¬ 
ments),  that  is  parallel  to  the  orientation  of  the  maximum  shear  stress. 
Consequently,  the  "soft"  elements  are  expected  to  slide  first;  and  if 
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then  the  "stiff"  elements  were  to  fail  by  a  tendency  to  separate,  the 
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above  observations  on  the  "failure”  sequence  would  be  confirmed.  Indeed, 
when  some  of  the  8  =  45°  elements  surrounding  the  8=0°  "stiff”  elements 
had  slid,  (due  to  roundoff  error  all  "soft”  elements  did  not  fail  at  the 
same  time,  but  were  sufficiently  close  to  sliding  to  be  considered  to 
fail  at  the  same  time),  negative  (tensile)  stresses  developed  in  all 
"stiff"  elements,  which  eventually  slid  until,  finally,  the  global  stiff¬ 
ness  matrix  became  singular  and  the  simulation  was  terminated.  The 
stress-strain  behavior  in  three  representative  elements  is  shown  in  Figs. 
59-61;  as  it  can  be  seen,  the  element  behavior  is  similar  to  the  one 
which  appears  in  Figs.  51-53.  Thus,  the  failure  of  the  aggregate  should 
be  attributed  to  this  "localization"  phenomenon,  which  is  similar  to  that 
observed  in  metals  (Dieter,  1976). 

This  "failure"  of  the  aggregate,  defined  here  by  the  sequence  of 
phenomena  previously  described  which  at  the  end  results  in  the  global 
stiffness  matrix  of  the  medium  becoming  singular,  is  associated  with  a 
generalized  tendency  of  the  particles  to  slide,  separate  and  rearrange 
themselves  into  more  stable  positions.  This  corresponds  roughly  to  the 
changes  in  geometry  occurring  in  actual  sands  at  the  threshold  strain, 

Yt  =  0.1  to  0.2  x  10”  3  (Dobry  et  al.,  1982),  as  verified  by  the  fact  that 
"failure"  of  the  medium  in  Fig.  50  occurs  at  a  shear  strain,  y  =  ej-e2  = 
0.16  x  10-3. 

In  addition  to  the  simulations  already  discussed  in  which  the  16  and 

64-element  media  were  sheared  at  constant  mean  stress  after  isotropically 

o  9 

compressing  them  with  aQ  =  1  Kg*/cm  ,  two  similar  runs  were  done  on 
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medium  2  (16-elements),  but  with  o0  =  0.5  Kg*/cm2  and  2  Kg*/cra2.  In 

these  two  runs,  the  medium  failed  at  tf  =  0.224  and  tf  *  0.896  Kg*/cm^, 

o 

respectively.  Therefore,  in  all  pure  shear  simulations  following  o0, 
o 

Tf/°o  53  0.448.  This  defines  the  “failure  surface"  shown  in  Fig.  54  for 

o 

the  aggregate,  where  the  corresponding  three  stress  paths  for  oQ  = 
constant  are  also  included.  The  observed,  and  proposed,  failure  surface 
is  a  cone  of  the  Von  Mises  type  with  the  slope  of  the  directrix  being 
0.448,  or  slightly  less  than  the  coefficient  of  intergranular  friction, 
f  *  0.5  (fig.  54). 


8.5  Discussion  of  Results 

As  it  can  be  seen  in  the  simulations  discussed  above,  the  media  used 
represent  well  some  aspects  of  the  behavior  of  actual,  uniform,  rounded 
sand,  by  being  isotropic  under  isotropic  pressure,  yielding  in  shear  and 
locking  under  hydrostatic  compression.  Therefore,  these  "random"  media 
have  the  desirable  properties  of  the  simple  cubic  array  without  its 
problematic  inherent  anisotropy  (see  Chapter  5). 

Since  the  elements  are  randomly  oriented,  discontinuities  in  the 
strain  field  are  expected  across  element  boundaries.  It  was  found  that 
the  maximum  strain  jump  at  the  nodes  was  one  to  two  orders  of  magnitude 
smaller  than  the  minimum  strain  value  in  the  16-element  media,  and  three 
orders  of  magnitude  smaller  in  the  64-element  medium.  This  is  thought  to 
be  satisfactory  in  both  cases  and  it  shows  that,  as  expected,  the  spatial 
accuracy  of  the  solution  improves  as  the  number  of  elements  increases. 

Even  though  the  macroscopic  behavior  in  Fig.  50  is  not  far  from 
being  linear  almost  until  "failure",  the  microscopic  behavior  is  not; 
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this  can  be  seen  in  Figs.  51-53  and  Appendix  E,  where  the  individual 


element  stress-strain  curves  are  shown  for  the  16-eleraent  medium  2.  The 


macroscopic  "linearity"  in  Fig.  50  is  thus  a  consequence  of  the  inter¬ 


action  between  all  the  elements. 


The  stress  field  inside  each  element  was  found  to  be  close  to 


uniform,  thus  verifying  a  basic  assumption  made  at  the  outset  of  this 


o  o  o  o 

study.  The  normalized  element  stresses  2o j j/ (01-02) »  2022/ (01-02)  and 


2oi2/(°l-02^ »  are  plotted  against  the  applied  shear  strain  ei~E2  for  all 


element  orientations  (8  orientations)  in  Figs.  G1-G8,  Appendix  G. 


Furthermore,  in  the  normalized  element  strain,  2eij/ei-e2>  is  plotted 
versus  the  applied  difference  in  principal  strain,  e[-£2  in  Fig®*  G9-G16 


of  the  same  Appendix  G.  It  can  be  observed  that  the  element  stress  and 


strain  fields  remain  almost  constant  during  loading. 


Moreover,  it  was  observed  that  the  stress  field  was  more  or  less 


uniform  throughout  the  medium  in  all  simulations,  except  very  close  to 


failure,  and  that  the  element  behavior  was  essentially  independent  of  its 


locations  and  a  function  only  of  the  orientation,  8,  of  that  element. 


(see  Fig.  49b  and  Figs.  E1-E17  in  Appendix  E) .  The  value  of  the  uniform 


stress,  o^j  for  each  element  was  approximately  the  applied  macroscopic 


stress,  o£j  resolved  in  that  particular  orientation,  ai j  “  a* j  ni,nj,  where 
o 

ojj  is  the  stress  applied  at  the  boundary  of  the  medium  and  n^nj  are 


direction  cosines.  This  was  true  in  all  runs.  This  is  illustrated  by 


Figs.  62-64,  which  show  the  variation  of  2aj  \/ (0 j-02) »  2a22^  ar|d 


2a 12/ (o i~a 2) t  versus  the  element  orientation,  8.  It  must  be  remembered 


that  0 1 1 ,  022  and  ai2  are  the  element  stresses  oriented  along  the  axes  of 
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the  cubic  array  (Fig.  45).  The  three  plots  in  Figs.  62-64  include  on, 
022  an<^  °12  calculated  for  all  elements,  (and  all  stress  increments  (with 
the  rectangles  indicating  the  range  of  values  for  each  increment),  in 
three  pure  shear  simulations  (a  =  0°,  45°,  and  90°),  and  the  curves  cor¬ 
respond  to  the  assumption  of  uniform  stress  field  throughout  the  medium. 
It  can  be  seen  that  the  assumption  of  the  independence  of  the  element 
stresses  on  element  location  is  fully  verified  by  the  plots.  It  should 
be  noted  that  this  is  identical  to  the  assumption  of  stress  uniformity  in 
the  early  work  of  Taylor  (1938),  and  very  similar  to  the  results  obtained 
by  Lin  (1964)  and  Budiansky  and  Wu  (1962)  for  an  assemblage  of 
elastoplastic  crystals. 


CHAPTER  9 


CYCLIC  LOADING  SIMULATIONS  OF  THE  AGGREGATE 

The  16-eleraent  media  were  subjected  to  a  number  of  stress  paths  with 
complete  stress  reversal  in  order  to  study  the  stress-strain  behavior  of 
the  aggregate  under  cyclic  loading,  and  also  to  compute  dynamic  proper¬ 
ties  such  as  the  damping  ratio. 

9.1  Cyclic  Isotropic  Compression 

The  (16-element)  medium  1  was  subjected  to  an  isotropic  stress  of 
o 

oQ  ■  3.0  Kg*/cmz  followed  by  a  complete  isotropic  cycle  with  amplitude  of 

,  9  0  , 

2.0  Kg*/cnr,  that  is  oQ  was  first  increased  to  5  Kg*/cm^  and  then 

decreased  to  1  Kg*/cra2,  and  then  back  to  5  KgVcm^.  The  stress-strain 

behavior  of  this  medium  is  shown  in  Fig.  65,  which  shows  a  nonlinear 

elastic  behavior  similar  to  that  of  actual  sand  (Ko  and  Scott,  1967). 

9.2  Isotropic  Compression  Followed  by  Cyclic  Shear  Loading 

o 

The  same  16-element  medium,  consolidated  isotropically  to  o0  =  1.0 
Kg*/cm2,  was  also  cycled  under  pure  shear  (a  =  45°)  conditions,  with  an 
amplitude  of  shear  stress,  tc  *  0.2,  0.3,  0.35,  0.4  and  0.43  Kg*/cm2. 

The  hysteresis  loops  for  xc  =  0.2,  tc  =  0.35,  xc  =  0.40  and  xc  3  0.43 
Kg*/cm^  appear  in  Figures  66a,  66b,  67a  and  67b  respectively. 

Computation  of  Dynamic  Properties  and  Wave  Velocities 

The  secant  shear  moduli,  G,  obtained  from  the  pure  shear  simulations 
at  different  a’s  for  the  constant  mean  stress  case  (Fig.  50),  were 
normalized  with  respect  to  the  shear  moduli  at  very  small  strains,  Gmax, 
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obtained  in  the  same  simulation,  and  the  corresponding  values  of  G/Gmax 
versus  shear  strain  are  plotted  in  Fig.  68a,  where  they  are  compared  to 
the  bounds  proposed  by  Seed  and  Idriss  (1970)  from  actual  tests  on  sands. 
The  corresponding  damping  ratio,  X,  was  obtained  from  the  loops  under 
cyclic  pure  shear  (Figs.  66a,  66b,  67 c,  67d)  with  the  expression: 

i-  AW 

2n  xy 

where  AW  is  the  area  of  the  stress-strain  loop,  x  and  y  are  the  maximum 
values  of  the  shear  stress  and  shear  strain,  respectively,  during  the 
cycle  (see  Fig.  69).  The  values  of  damping  ratio  obtained  in  cyclic  pure 
shear  are  plotted  agains  cyclic  shear  strain  y  in  Fig.  68b.  Figure  68 
corresponds  to  simulations  done  at  a  mean  stress  of  1.0  Kg*/cm^.  Figure 
68  also  includes  the  corresponding  curves  for  G/Gmax  and  damping  ratio 
for  the  case  of  pure  shear  along  the  slip  planes  of  a  simple  cubic  array, 
computed  by  Dobry  et  al.  (1982).  Since  the  slope  of  the  G/Graax  and 
damping  ratio  versus  y  curves  in  Fig.  68  for  the  aggregate  studied  herein 
is  flatter  than  that  of  the  simple  cubic  array,  its  behavior  is  more 
realistic  when  compared  to  actual  sand.  This  "stiffer"  behavior  of  the 
aggregate,  as  compared  to  that  of  a  single  cubic  array,  should  be 
attributed  to  the  interaction  between  the  elements  of  the  medium. 

9.3  Constrained  Moduli  of  the  Aggregate  and  P-Wave 
Velocities  Under  Biaxial  Compression 

The  small  strain  constrained  moduli,  j ,  were  obtained  from  the 


"triaxial”  simulation  previously  discussed  (Isotropic  Compression 
Followed  by  Biaxial  Compression,  Section  8.3).  The  values  of  Djj  were 
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computed  in  both  principal  (vertical  and  horizontal)  directions.  In 

order  to  increase  the  accuracy  of  the  calculation  at  these  small  stress 

and  strain  increments,  and  to  fully  account  for  the  stress-induced 

anisotropy  of  the  aggregate,  the  64-element  medium  (medium  4)  was  used. 

o  o 

The  medium  was  loaded  to  the  desired  biaxial  stress  ratio,  K  =  01/02 >  ar|d 

then  very  small  stress  increments  with  appropriate  signs  were  applied, 
o  o 

AO|,  Ao2»  the  corresponding  differences  in  strain  were  computed  and, 

finally,  the  small  strain  constrained  moduli  of  the  medium,  Du,  D22  were 

computed  in  both  directions.  The  results  of  this  simulation  are  shown  in 

Figs.  70a  and  70b  as  plots  of  the  normalized  constrained  moduli  of  the 
(K)  (1)  (K)  (1)  0  0 

medium,  D22  /022  ar>d  D1 1  /°11  versus  the  stress  ratio  K  *  01/02*  where 

00  (1) 

Dj£  is  the  constrained  modulus  at  a  given  K,  and  Du  the  corresponding 

o 

constrained  modulus  under  the  initial  isotropic  stress,  oQ.  The  same 
figure  also  includes  data  points  from  a  number  of  measurements  on  sand 
obtained  in  the  large  cubic  facility  at  the  University  of  Texas  at  Austin 
(Kopperman  et  al. ,  1982),  which  were  performed  during  a  test  with 
conditions  similar  to  those  assumed  in  this  numerical  simulation.  The 
agreement  between  experimental  results  and  numerical  simulations  in  fig. 
70a  and  70b  is  excellent.  Consequently,  the  main  conclusion  obtained 
from  the  University  of  Texas  laboratory  results,  that  the  P-wave  velocity 
propagating  along  a  principal  stress  direction  is  only  a  function  of  the 
value  of  that  principal  stress  is  fully  predicted  by  the  numerical  experi¬ 
ment.  Therefore,  as  previously  hypothesized  by  the  authors,  this  effect 
is  due  to  the  particulate  nature  of  the  soil,  and  its  analytical  modeling 
necessarily  requires  taking  into  account  this  particulate  nature. 


CHAPTER  10 


CONCLUSION 

A  particulate  mechanics  model  has  been  developed  for  describing  the 

elastic  response  of  an  assemblage  of  identical,  elastic,  rough  spheres  of 

arbitrary  macroscopic  porosity,  n,  subjected  to  an  isotropic  pressure, 
o 

oQ.  The  model  is  based  on  the  Mindlin-Deresiewicz  theory  of  bodies  in 
contact  and  takes  into  account  the  spatial  variation  of  porosity.  The 
model  assumes  that  the  assemblage  is  composed  of  random  clusters  of 
several  regular  arrays  of  various  porosities  and  it  computes  the  macro¬ 
scopic  moduli  by  means  of  the  Self  Consistent  Method. 

The  predictions  of  the  model,  specialized  for  the  case  of  quartz 
spheres,  were  compared  to  measurements  of  the  shear  modulus,  Gmax,  on 
uniform  quartz  sands,  with  good  qualitative  agreement;  however,  the 
analytical  "sands"  were  as  much  as  3.5  times  stiffer  than  the  actual 
soils.  This  is  explained  by  the  fact  that  sands  have  less  effective 
contacts  per  grain  than  theoretically  predicted  for  a  given  porosity,  n. 
However,  when  the  sand  is  prestrained  by  millions  of  shearing  cycles 
slightly  above  the  threshold  strain,  the  shear  modulus  approaches  the 
theoretical  value  without  changing  the  macroscopic  porosity,  n,  as  the 
theoretical  number  of  contacts  is  slowly  realized.  Thus,  the  model 
exhibits  excellent  agreement  with  results  on  heavily  prestrained  sand, 
and  it  also  provides  upper  bounds  for  small  strain  shear  and  bulk  moduli 
for  rounded,  uniform  sands. 

Finally,  a  nonlinear  finite  element  formulation  is  established  which 
interprets  the  behavior  of  granular  soil  at  small  strains.  This 
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formulation  also  considers  the  soil  as  an  aggregate  of  different  packings 
of  equal  spheres,  but  this  time,  for  simplicity,  only  the  simple  cubic 
packing  is  used.  The  inherent  anisotropy  of  the  simple  cubic  packing  is 
eliminated  by  orienting  the  same  packing  randomly  in  several  directions. 
The  model  incorporates  the  nonlinear  constitutive  model  CONTACT  for 
describing  the  force-deformation  behavior  at  the  contact,  and  results  are 
obtained  for  several  stress  paths.  Specifically,  the  recent  compressional 
wave  velocity  measurements  on  anisotropically  consolidated,  dry  sand  are 
simulated.  These  experimental  measurements  have  shown  that  the  P-wave 
velocity  depends  only  on  the  principal  stress  on  the  direction  of  wave 
propagation;  this  was  predicted  by  the  numerical  simulation.  Further¬ 
more,  the  comparison  of  the  normalized  constrained  moduli  between  the 
finite  element  model  and  the  measurements  is  excellent,  thus  validating 
the  micromechanics  approach  used  to  the  problem. 

Both  phenomena  described  above,  the  increase  in  the  stiffness  of  a 
granular  soil  and  the  dependence  of  the  P-wave  velocity  only  on  the 
principal  stress  in  the  direction  of  wave  propagation,  had  not  previously 
been  fully  interpreted  or  modelled  analytically.  Their  interpretation 
and  modeling  here  was  made  possible  only  through  a  micromechanical 
approach,  which  takes  into  account  the  particulate  nature  of  the  granular 
soil  responsible  for  the  above  phenomena.  Therefore,  the  method  of 
modeling  rounded,  uniform  granular  soil  as  an  aggregate  of  regular 
packings  of  equal,  elastic,  rough  spheres  leads  to  simple  and  accurate 


results 
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(1.1.21 

0. )985 

0.6625 

21 

12.2.4) 

0. 1954 

0.6540  ~ 1 

22 

1  3, 2. 3) [2. 4. 21:11. 6. 11 

0.  1954 

0-6529 

2) 

[4,0.41 

0. 3198 

0.4702 

24 

9 

[l.*.4|  . 

0.  1866 

0.6)0) 

25 

[2,1,21=11,6.21 

0  3520 

0.54)2 

26 

10 

12.4.41=11.6. li 

0.  3)4) 

0.5022 

27 

[2,5.11 

0. 3127 

0.4550 

28 

14. 2. 41=12.6.21 

0.  3019 

0  4)27 

29 

11 

(2.6. JJ 

0281) 

0  )908 

30 

12 

( 3.6,  3) 7 [4,4 ,4 ) 

0.2595 

0 . 3504 

_ 

31 

Table  2.  Feasible  Regular  Arrays  or  "Cells"  (Shahinpoor, 
1981).  Coordination  No.  ■  N  -  No.  of  Contacts 
per  Sphere.  [u,m,ij  gives  No.  of  Contacts  of 
Spheres  with  Layer  Above,  Same  Layer  and  Layer 
Below:  N  ■  u  +  m  +  i. 


Young's  Modulus 

Es  -  11  x  106  psl 

Poisson's  Ratio 

vs  -  0.15 

Coefficient  of 
Intergranular  Friction 

f  -  0.5 

Table  3  Properties  of  Quartz  Used  in  this  Report 


{Note  that  Lambe  and  Whitman  (1969)  reported  a  different  value 

■  0.31  for  quartz,  which  in  turn  was  used  for  the  calculations 
in  Dobry  et  al.  (1982).  The  lower  value  vs  *  0.15  used  herein 
is  more  realistic  and  was  obtained  from  White  (1964)  and  Ko  and 
Scott  (1967)] 


Shear  Wave  Velocity,  V  ,  mps 
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Variation  of  Vg  with  Principal  Stress  in  (a)  Direction 
of  Wave  Propagation,  (b)  Direction  of  Polarization 
(Particle  Motion)  and  (c)  Out-of-Plane  Direction 
(Roesler,  1979) 


Figure  1. 
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Stress  Perpendicular  to  Wave  Propagation,  o  ,  psi 
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Figure  2.  Effect  on  Vp  of  Principal  Stress,  ac,  Perpendicular  to 
Wave  Propagation  for  Biaxial  Loading  Case  (Kopperman 
et  al. ,  1982) 
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Figure  3.  Stress-Strain  Curves  for  Monotonic  Loading  of  Dry 
Granular  Soils. 
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Figure  4.  Stress-Strain  Hysteresis  Loops  for  Reversed  Loading 
(Hardin  and  Drnevich,  1972)  . 
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Figure  9.  Elastic  Spheres  Under  Normal  and  Tangential  Loads 


Figure  10.  Tangential  Force-Displacement  Relation  for  N 
Constant,  T  Increasing  (Dobry  et  al,  1982) 
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(a)  simple  cubic  (b)cubical  tetrahedral 


(c)  tetragonal  (d)  pyramidal  (a)  tetrahedral 
sphenoidal 


Figure  11.  Regular  Arrays  of  Equal  Spheres  (Deresiewicz ,  1958) 
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Figure  16.  Force-Deformation  Behavior  of  Two  Elastic  Rough 
Spheres  in  Contact:  Analytical  Solutions  for 
(a)  Normal  Force  Increasing,  Tangential  Force 
Increasing;  (b)  Normal  Force  Decreasing, 
Tangential  Force  Increasing;  (c)  Normal  Force 
Increasing,  Tangential  Force  Decreasing; 

(d)  Normal  Force  Decreasing,  Tangential  Force 
Decreasing  (after  Mindlin  and  Deresiewicz,  1953) 


TANGENTIAL  DISPLACEMENT,  <5(CM)(x  10~3) 


Numerical  Simulation  of  the  Force-Deformation  Behavior 
for  Two  Equal,  Elastic,  Rough  Spheres  in  Contact: 
Load-Displacement  Relation  for  an  Oscillating  Oblique 
Force  with  dT/dN  <  f 


Figure  18 
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'igure  20. 


Numerical  Simulation  of  the  Force-Deformation  Behavior 
for  Two  Equal,  Elastic,  Rough  Spheres  in  Contact: 
Load-Displacement  Relation  for  the  Case  in  Which  the 
Normal  Force  Decreases  and  the  Tangential  Force 
Increases 
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TANGENTIAL  DISPLACEMENT,  <5(CM)(x  10~3) 


Figure  22.  Numerical  Simulation  of  the  Force-Deformation  Behavior 
for  Two  Equal,  Elastic,  Rough  Spheres  in  Contact: 
Load-Displacement  Relation  for  the  Case  in  which  the 
Normal  Force  Decreases  and  the  Tangential  Force 
Decreases 
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a)  Body-centered  cubic  unit  cell  b)  Body-centered  cubic  structure 


Figure  25.  (a)  Body-Centered  Cubic  Array  of  Equal  Spheres; 

(b)  Body-Centered  Cubic  Array  (after  Van  Vlack, 
1964). 
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a)  Face-centered  cubic  unit  cell  b)  Face-centered  cubic  structure 


Figure  26.  (a)  Face-Centered  Cubic  Array  of  Equal  Spheres; 

(b)  Face-Centered  Cubic  Array  (after  Van  Vlack, 
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Figure  29.  Secant  Shear  Modulus  Versus  Shear  Strain:  Comparison 
Between  Calculated  G/Gmax  for  the  Simple  Cubic  Array 
and  Experimental  Range  for  Sand  (Dobry  et  al., 

1982) 
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Figure  30. 
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Influence  of  Stress  Ratio,  K  -  a]/a3  on  GLax: 
Comparison  Between  Prediction  From  Simple  C^bic  Array 
and  Experimental  Range  for  Sand  X 
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Figure  32. 


Poisson's  Ratio  of  Regular  Cubic  Arrays,  v,  as  a 
Function  of  the  Poisson's  Ratio  of  the  Spheres,  vs 
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Triaxial  Loading:  Obliquity,  aC  Fail1 

Versus  the  Coefficient  of  Intergranular  Fric 
for  Three  Regular  Cubic  Arrays. 


Figure  33 
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Figure  34.  Shear  Wave  Velocity,  Vs,  Versus  Stress  Ratio,  K, 
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Figure  46.  16-Element  Media  Used  in  the  Simulation  and  Angle 
of  Orientation  of  Each  Element. 
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Figure  47.  Stress-Strain  Behavior,  (c),  of  a  Simple  Cubic  Array 
of  Quartz  Spheres  Subjected  to  Biaxial  Stress, 
0U*1.61,  O22“2.0  KgVcra^,  Followed  by  Pure  Shear, 
oi2»  to  Failure,  (a)  and  (b)  Portray  the  Correspond¬ 
ing  Force-Deformation  Behavior  of  the  "Weak"  and 
"Strong”  Contact,  Respectively. 
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Figure  48.  Medium  1:  Stress-Strain  Behavior  for  Hydrostatic 
Compression  Loading  starting  from  <jo=3.0  Kg*/cra2. 


(a)  Orientation  of  the  Applied  Principal  Stresses 
During  Pure  Shear  Loading  with  Mean  Stress  Constant 
ao=0.56(a°  +  o^^Constant;  (b)  Orientation  of  an 
Element 
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Figure  50. 


o  o 

Difference  in  Applied  Principal  Stresses,  01-02, 
Versus  Shear  Strain,  Y=ei~e2»  Calculated  for  all 
Media  (16  nd  64  elements),  and  for  Values  of  a=0°, 
22.5°,  45°,  67.5°,  and  90°  (MX=medium  X, 
AYY-inclination  of  01  is  YY  degrees). 
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Figure  55. 
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Medium  l:  Stress-Strain  Curve  for  the  Aggregate  for 
Compression  Loading  with  Variable  Mean  Stress 


Following  Isotropic  Pressure,  Oq*1  Kg*/cm^. 
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Figure  57.  Medium  1:  Compression  with  Variable  Mean  Stress 

(a£*1.0  Kg*/cnr)  Stress-Strain  Behavior  of  Element 
E2  Oriented  at  8=50°. 
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Figure  58.  Poisson’s  Ratio  for  the  Aggregate  Calculated  from 
the  Compression  with  Variable  Mean  Stress 
Simulation. 
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Figure  62.  Medium  2:  Compression  with  Constant  Mean  Stress 


(<7o*1.0  Kg*/cnr).  Normalized  Element  Normal  Stress, 
ojl/(o°-a2)  Versus  Angle  of  Orientation  B,  of  the 
Element  for  Three  Inclinations  of  Principal  Stress, 
a«0,  45  and  90  degrees. 
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Figure  63.  Medium  2:  Compression  with  Constant  Mean  Stress 

(oq-I.O  K£*/cm^).  Normalized  Element  Normal  Stress, 
a22/(<Jl~ol)  Versus  Angle  of  Orientation  8»  of  the 
Element  for  Three  Inclinations  of  Principal  Stress, 
a-0,  45  and  90  degrees. 
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Figure  64.  Medium  2:  Compression  with  Constant  Mean  Stress 

(ao*1.0  i^*/cm2).  Normalized  Element  Normal  Stress, 
0l2/(ai~02)  Versus  Angle  of  Orientation  8,  of  the 
Element  for  Three  Inclinations  of  Principal  Stress, 
a«0,  45  and  90  degrees. 
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Figure  68.  (a)  Normalized  Shear  Modulus,  G/Graax;  (b)  Damping 

Ratio,  X;  both  Versus  Shear  Strain,  y.  Comparison 
Between  Values  Measured  in  Sand,  and  Those  Calculated 
with  the  Simulated  Aggregate  and  the  Simple  Cubic 
Array  in  Pure  Shear  . 


O  Experiment  (Kopperman  et  al  1982) 
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APPENDIX  A 

STRESS-STRAIN  RELATIONS  FOR  A  BODY  CENTERED  CUBIC  ARRAY 

Following  a  procedure  similar  to  that  used  for  the  simple  cubic  array 
by  the  authors  and  by  others,  incremental  stress-strain  relations  for 
various  states  of  initial  stress  can  be  obtained  for  the  body  centered 
cubic  array. 


A.l  Relation  Between  Stress  and  Contact  Forces 
Consider  first  a  medium  composed  of  identical  spheres,  Fig.  25b, 
arranged  in  a  body  centered  cubic  array.  Take  as  an  element  of  the  medium 
the  cube  shown  in  Fig.  25c.  This  "elementary"  cube  (or  representative 
volume).  Fig.  AI,  is  chosen  to  contain  a  sufficient  portion  of  the 
medium  to  define  the  arrangement.  Clearly,  each  sphere  in  the  medium  is 
in  contact  with  8  other  spheres. 

Increments  of  the  force  dPjj  act  on  the  forces  of  the  cube,  Fig.  Al , 
and  they  are  assumed  to  be  distributed  among  the  spheres  in  proportion  to 
their  stiffness,  that  is  to  their  section  exposed  on  the  faces  of  the 
cube. 


The  incremental  stress  components  are  defined  as  follows: 
d°ij  “^dPij/(-Y-) 


(Al)  j 


where  iSSl 
3 


is  the  gross  area  of  the  face  of  the  cube. 


At  each  contact  between  spheres,  the  normal  and  tangential  components 
of  the  incremental  force  are  again  designated  by  dNjj ;  with  dNj^  being  the 
normal  component  and  dNjj  being  two  the  tangential  components. 


Once  more,  the  first  step  in  deriving  the  incremental  stress-strain 
relationships  is  to  define  the  expressions  for  the  increments  in  the 
forces  at  the  contacts  between  the  spheres  in  the  cube.  Fig.  25a,  dNjj,  in 
terms  of  the  increments  of  the  applied  stress  dajj.  However,  since  this 
array  is  statically  determinate  for  initial  isotropic  and  transversely 
isotropic  triaxial  loading,  only  the  equilibrium  conditions  are  sufficient 
for  the  solution  of  this  subproblem.  However,  this  case  is  much  more 
involved  than  the  simple  cubic  array  and  tedious  calculations  have  to  be 
performed. 

Fig.  A2  shows  one  octant  of  a  sphere  at  the  apex  H  as  well  as  the 
point  of  contact  with  the  "central"  sphere  and  the  applied  and  contact 
forces.  This  octant  of  the  sphere  at  H  will  be  treated  as  the  representa¬ 
tive  octant.  Now  equilibrium  equations  have  to  be  written  for  each 
spherical  octant  at  every  apex  and  the  contact  forces  will  have  to  be 
solved  for  each  case  separately.  For  example  for  apex  H,  Fig.  A2,  the 
equilibrium  equations  are:* 


*  The  primed  symbols  refer  to  the  local  coordinate  system. 
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I  Fx,  a  0  =>  —  ^33 - —  dN^ - —  dN^2  +  j  dP13  +  j  dP33  +  j  dP23  = 


the  solution  of  which  is: 


dN31 - -  (dPn  -  dP33  +  dP12  -  dP23) 

O 


dN^  =  (dPn  -  2dP22  +  dP33  -  dP12  +  2dP13  -  dP23)  (A6) 

24 


dNss  =  —  —  (dP^i  +  dP22  +  dP33  +  2dP^2  +  2dP}3  +  2dP23)  (A7) 


At  this  point  the  state  of  stress  can  be  defined  and  the  constitutive 
law  may  be  determined  for  each  case  (isotropic  or  cross-anisotropic 
triaxial  confinement). 


A. 2  Isotropic  State  of  Stress: 


Applying  increments  of  force  along  the  three  principal  directions. 


dPn,  dP22i  dP33  and  one  at  a  time,  on  top  of  the  isotropic  confining 
stress,  an  inremental  force-deformation  relationship  can  be  developed. 
For  example,  in  the  case  of  application  of  dP^j  (Fig.  Al)  we  have: 


"  (i  ^  +  i Ct)  dPu 


d«22  -  +-  <cn  -  Ct>  dPll 


d«33  *  ~  (cn  -  Ct>  dPll 
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where  Cjj,  Ct  are  the  normal  and  tangential  Compliances  at  the  contact. 
Similarity  applying  dP22: 

(All) 


d«ll  -  Cn  Ct)  dP22 


d«22  =  "  (7  °n  +  J  Ct)  dP22 


dd3 3  -  C„  Ct)  dP22 


(A12) 


(A13) 


To  determine  now  the  relation  between  changes  in  angle  and  forces  we 
have  to  look  at  the  difference  between  displacements  (Fig.  Al).  For 
example,  dPj2  being  applied,  we  have: 


dY23 


*»|a  -  633U  ,  622U  -  822!d 


abcc 


abcc 


etc. 


(A14) 


where  af,cc  is  the  length  of  edge  of  the  cube.  Now: 


633|a  “  633|Ie  =  2(~  533  +  ~  6 13  +  623) 


622 


-  622  =  0  etc. 

A  D 


(A15) 


(A16) 


Consequently : 

dYl  2  -  (|  ^  +1  Ct)  dP12  (A17) 

Evaluating  the  compliances 

l-vs 

Cn  “  Tn -  which  yields 

n  2Gsa0 


n 
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(l-vs)2/3  ^  2 

(4VTGs20o)l/3  R 


(A18) 


Similarily 


2(l-vs)l/3  (4VTgs20q)1/3 


(A19) 


in  the  case  that  vg  *  0 


d«ll  =  C„  +  j  Ct)  dPH 


Substituting  for  Cn,  Ct  eqns.  (A18,  A19)  we  obtain: 


6  (  4V~Tgs  2  o0 )  1  /  3  « 


-  dPn  +  - 

D  1  *  £ 


~ — -  - - - -  dPn 

6  (l-vs)l/3  (4V^TGs2ao)l/3  R 

(A20) 


then: 


deu  =  — —  — - —  [0-vs)2/3  + - ~7:]  d°ll 

373  (4V3Gg2ao)l/3  (l-vs)1/3 


(A21) 


3v/3  (4VTg|  CTo)i/3 

dan- - -  d£ll 

2[(l-vs)2/3  + - ^-r] 

(i-vs)1/3 


(A22) 


Similarily 


de22  -  - - - -  [2(  l~vs) 2/3 - — - ]  dou 

3Vr(4V^Gg2ao)l/3  CI-vs)l/3 


(A23) 
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d£33  =  ds22 

Also  from  eqns.  (A17,  A18)  and  ( A 1 9 )  we  have  that 


(A24) 


de  1  2 - “ 


3  VT  (4VTgs20o)1/3 


[(l-vs)2/3  +- - — — ]  doi2 

4  ( l“vs) 


(A25) 


Finally,  the  Incremental  stress-strain  law  for  the  isotropic  case  may  be 
written  as  follows 

dell  =  cllll  doll  +  C1 122  do22  +  c 1 l 33  da33 

de 22  *  C1I22  d^ll  +  c2222  d^22  +  c2233  d°33 

de33  =  Cj^33  dan  +  ^3322  dar22  +  C3333  da33  etc.  (A26) 

where  are  the  expressions  (compliances)  in  the  stress-strain 


le  1 1 


l£  22 


d< 
de33 
de  1 2 
de  1 3 
de23 


1 

! 


example 

,  in  eqn 

.  (A21), 

(A22). 

In  matrix 

form 

-i 

r  i 

« 

* 

cl  1 1 1 

C1 122 

C1 133 

0 

0 

0 

don 

1 

| 

c  2  2 1 1 

c2  2  2  2 

c2233 

0 

0 

0 

do22 

1 

c  3  3 1 1 

c3322 

c3  3  33 

0 

0 

0 

do33 

* 

* 

0 

0 

0 

C1 2 12 

0 

0 

dai  2 

0 

0 

0 

0 

C1  3  1  3 

0 

daj3 

J 

i 

0 

0 

0 

0 

0  C2323 

da2  3 

(A27) 

p 

> 

i 

j 

i 

i 

i 

bcc  array  to  be 

isotropic 

under 

isotropic 

loading 

the 

following  condition  must  be  satisfied: 


V' 
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C1 212  =  c1313  =  c2323  *  CU11  -  C1 1 22  (A28> 

c2 2 1 1  3  c  1 1 22  =  C1L33  =  c2233  =  c33 1 1  =  c3322  =  CU22  (A29> 

The  above  conditions  are  satisfied  only  when  vs  =  0;  furthermore,  in  this 

case  of  the  bcc  array,  as  in  the  sc  array,  C j  ^ 22  =  0  an8  the  compliance 
matrix  is  diagonal. 

In  this  case: 

2  1  1/3 

C1 1 1 1  3  c2222  3  c3333  3  -  ( - - - )  (A30> 

V3“  ^VT"gs2  o0 


Since  the  compliance  matrix  is  diagonal,  its  inverse,  the  stiffness  matrix 
is  easily  computed  by  inverting  each  term:  i.e. 


(A31) 


and  clearly  the  shear  modulus  of  the  array, 


G,  is  G  (4V3-  Gs2 


(A32) 


At  this  point  the  Poisson's  Ratio  of  the  bcc  array  may  be  computed  as 
follows: 


e33 

1  1  e22 1 

ei  1 

r  I £  1 1 

(l- 

,2/3  _ 

2  — Vg 

1/3 

vs> 

2(1-Vg) 

(1- 

,2/3  X 
vs)  + 

2-Vg 

(l-vs) 

1/3 

(A33) 


i 

I 

i 


(A34) 


For  different  values  of  vg ,  the  Poisson's  Ratio  of  the  bcc  array,  vt>cc> 
may  be  computed;  a  plot  of  v^cc  versus  the  Poisson's  ratio  of  the  spheres, 
vs,  appears  in  Fig.  (32)  together  with  the  variation  of  vsc  and  vfcc  with 
vs  for  easy  comparison. 


A. 3  Transversely  Isotropic  State  of  Stress  (Triaxial  Loading) 

As  in  the  case  of  the  simple  Cubic  Array,  the  application  of  an 
anisotropic  stress  increment  will  result  in  a  variation  of  the  contact 
forces  and,  consequently,  of  the  corresponding  contact  stiffnesses. 
Therefore,  in  order  to  obtain  the  stress  strain  relationships,  the 
derivation  must  be  performed  once  more  distinguishing  between  compliances 
at  contacts  with  different  loading  histories.* 


*The  computation  of  the  compliances  for  the  case  of  anisotropic  loading 
for  both  the  sc  and  the  bcc  arrays  has  been  done  in  order  for  the  results 
to  be  used  only  for  the  special  case  of  wave  propagation.  This  way,  the 
load  was  assumed  to  reverse  direction,  and  for  this  the  elastic  tangen¬ 
tial  compliances  were  used.  In  the  general  case,  the  load  could  either 
increase  or  decrease  monotonical Ly  and  different  compliances  in  each  case 
would  apply. 
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Consider  now  the  case  of  cross-anisotropic  stress  ("triaxial  test") 
imposed  after  the  array  has  been  subjected  to  an  initial  hydrostatic 
stress.  That  is  (Fig.  A3): 


at  t  -  tQ:  ou  =  022  =  033  =  o0 


t  =  tL:  on  =  033  -  aQ 


a22  =  ao  + 


(A35) 

(A36) 

(A37) 


The  contact  forces  during  this  Anisotropic  Loading  are: 
dN31  =  ±  ~ 8~  ^Pl1  ~  d?33^ 


(A38) 


•  VS 

dN32  =  *“24  fdPl1  ~  2d?2  +  dP33l 


(A39) 


.  VT" 

dN33  =  i  ~i2  [dPll  ~  dP22  +  dp33] 


(A40) 


in  the  case  of  transversely  isotropic  loading  the  above  equations  simplify 
to 


dN3 1=0 


j  n •  VT 

dN32  =  ±  -  dPa 

12 


1  VT" 

dN33  =  ±  -  dPa 

12 


(A41) 


(A42) 


(A43) 


This  is  the  case  of  the  two  spheres  in  contact  where  the  normal  force  is 
increasing  from  N0  to  N0  +  N*  and  the  tangential  force  from  0  to  T*,  with 
d  T 

6  • —  >  f  (Mindlin  and  Deresiewicz,  1953) 
dN 


■t.n 


>3 

$ 

'/> 

I 


'.N'J 


(usually  0.5  _<  f  _<  0.8  for  sands) 

in  this  case  the  Normal  Compliance,  Cn,  is: 

l-vs 


O  3( 1-Vg) 

where  a^  =  — -  R(N0  +  N) 


and  finally 


.  3(l-vs)  „  i  ,aas! 

“3  «o  [1  *7  t)l 


2G, 


s 


then 


(^L-]l/3 

L4VTg  2. 


-1/3 

i 


The  Tangential  Compliance,  Ct,  is: 


2-v  L*  L  ^  ^ 

ct  -  4^f  t-9  *a.e)(l-(i-0)  } 


.  _  f  .  T  ,  T*  T* 

where  9  =  — ,  L  =  -  and  L  =  - 

8  fNo  fNo 


However,  in  this  case,  since  the  load  decrement  is  small,  then  L 


(A45) 


(A46) 


(A47) 


(A48) 


*  T 

-  L  ♦ 


therefore 


177 


2-vs 
Ct  "  4G^ 


in  this  case 


(A49) 


Ct  =— ^ -  [ - i - ]’/3  \l  [^)]  1/3~ 

2(l-vs)1/3  4 VT~Gs2oq  3  °°  R 


(A50) 


The  constitutive  law  for  this  particular  loading  may  be  developed  in  the 


same  manner  as  in  the  case  of  the  isotropic  loading,  i.e. 


d«22  *  Cn  +  ~  ct]  dP22 


f inally 


(A51) 


2  ,  1  ^/3r,  .  1  1/3rn  ,2/3  2“vs 

-  [1+-(~)J  [0-vs)  + - — ]  do22 

3VT~  4VTlls2a0  3  °  (l-vs)1/3 


(A52) 


for  vs  ■  0. 


•p  -  I/3  .  0  -1/3 

de22  = -  (  — )  [1+T^r)l  da22 

vr  ^vtgs20o  3  °° 


(A33) 


in  the  case  of  isotropic  loading  and  vs  *  0,  the  above  reduces  to 


2  ,  1 

de22  ”  -  ( - 1  do22 

VT"  ^~GS2  0o 


(A54) 


which  is  the  same  obtained  before,  eqn.  (A30).  Analogously 
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1  1/3  1  0a  -1/3 

(-  ■.  )  [1+i(^)l  [<‘-vsr 

U\r^a2„  3  ao 


VT  4VT"gs20o' 


2( l-vs) 1/3 


]  d°22 


(A55) 


For  vs  =  0  this  reduces  to  zero. 


The  Shear  Compliance  is  computed  as  follows: 


d«ij  =  [f  °n  Ct]  dPij 


(A56) 


y  1  */3  1  0_  -1/3  9/0  2-Vc 

—  ( - 5 - )  [i+4(— )]  [4<i-vs)2/3  * — - ]  d0 

3V5-  4VT"G32o0  3  °  2<l-.Js)l/3  J 


(A57) 


For  us  =»  0  this  reduces  to 


2  1  1/3  1  Oa  _1/3 

dc  ,-L.{  1  )  day 

wr  *vr gs20o  3  a° 


(A58) 


For  aa  »  0  this  reduces  to  eqn.  (A31).  Therefore,  for  the  case  of  vs  H 


2  1  1/3  la  _1/3 

C1 1 1 1  *  c2222  =  c3333  = -  ( - 7 — )  \ 1  +  ~  (t2-)  ]  * 

3VT  ^VTgs2CTo  3  o0 


*  r,,  ,2/3  2-vs 

[(1-Vg)  +- - —  ] 


(l-v,,)1/3' 


(A59) 
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-2  1  1/3  1  oa  "1/3 

C1 122  *  c2233  "  c33 1 1  - - - — )  [*  +~  (— )  ]  * 

vr  «vr  cs2ao  3  o0 


[(l-vs) 


2/3  2~vs 


2(l-vs)1/3 


(A60) 


1  1  1/3  1  Oa  "1/3 

G 1 2 1 2  =  c1313  =  c2323  Z - - — )  [*  +~  (~ )  ]  * 

3 VT  bVT  Gs2ao  3  °0 


*  r /  / .  ,2/3  .  2-vs 

*  [4<l-vs)  +- - — ] 

(l-vs)l/3 


(A61) 


"1 


In  this  case  the  material  will  be  isotropic  under  cross  anisotropic  load¬ 
ing  only  if  in  the  compliance  matrix 


C11U  ~  122  =*  C1212 


Performing  the  calculations  we  see  that  indeed. 


i  i  I/3  ,  a  -1/3 

:lill_c1122=  - {— - - — )  [1+“(~)]  [40-vs)  + 

3vT4VT  Gs2ao  3  °o 


2-Vs 


(l-v.)‘/3 


]  “c  1 2 1 2 


Therefore,  the  bcc  array  is  isotropic  under  cross  anisotropic  loading*; 
this  is  a  serious  deficiency  of  the  model  and  should  be  attributed  to  the 
symmetry  of  the  array. 

In  the  case  that  vs  =  0,  the  array  is  still  isotropic  but  this  time 
the  compliance  matrix  is  again  diagonal: 


*f 


or  the  conditions  specified  previously. 
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Simulation  of  Triaxial  and  Pure  Shear  Loading 
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APPENDIX  B 


SIMULATION  OF  TRI AXIAL  AND  PURE  SHEAR  LOADING  IN  CUBIC  ARRAYS 


The  three  regular  cubic  arrays  will  be  subjected  to  a  finite  loading 


in  such  a  way  so  as  to  cause  failure  when  applied  along  a  principal 


direction.  In  the  case  of  the  simple  cubic  array,  this  corresponds  to  a 


pure  shear  loading,  since  for  a  triaxial  loading  the  array  locks. 


Conversely,  the  body  centered  cubic  and  the  face  centered  cubic  array  will 


be  subjected  only  to  a  triaxial  loading,  since  for  a  pure  shear  loading 


the  array  also  lock. 


In  order  to  obtain  a  finite  stress  strain  relation  in  every  case,  the 


compliances  need  to  be  integrated  along  the  loading  path.  Once  this  is 


done,  finite  displacements  are  computed  and  then  the  strains  are  obtained 


in  the  same  manner  as  in  the  infinitesimal  constitutive  laws. 


Finally,  the  load  is  increased  until  failure  (gross  sliding)  occurs 


in  the  array;  in  the  statically  determinate  arrays  gross  sliding  at  a 


contact  translates  into  failure  of  the  cubic  array.  In  the  statically 


indeterminate  face  centered  cubic  (fee)  array,  the  medium  does  not  fail 


immediately,  but  first  the  number  of  contacts  reduces  from  12  to  8  when 


the  fee  array  becomes  statically  determinate,  and  then  sliding  at  one 


contact  becomes  failure. 


B.l  The  Simple  Cubic  Array  Subiected  to  Pure  Shear  Loadin 


Consider  the  Simple  Cubic  Array  shown  in  Fig.  17  and  consider  a 


force  T  acting  in  the  X[  direction  (Fig.  Al).  The  s.c.  array  is  subjected 


to  an  initial  isotropic  force  Nc  and  the  value  of  T  increases  monotonical ly 


from  zero  to  T*,  where  T*  is  the  value  of  T  that  causes  failure  in  the 
array,  while  NQ  remains  constant.  In  this  case  the  tangential  displace¬ 
ment  is  given  by  Mindlin  (1949): 


6 


3(2-us) 
8Gg  a 


f  N0 


[1  -  (1 


JL  ] 

f  N_ 


2/3 


(Bl) 


where  a 


3 


3(l-vs)N0R 

8GS 


is  the  radius  of  contact. 

Now  y  =■  6/R;  substituting  the  expression  for  a^  into  eqn.  (Bl)  and  after 
transforming  the  forces  into  stresses  we  obtain: 


Y 


2  1/3 

3<2~vs)  *  f  r  qo  1 

(l-vs)l/3  12Gs2 


[1  -  (1 


T 

f  00 


2/3 

)  ] 


(B2) 


The  above  equation  has  been  plotted  in  Fig.  28  for  different  values  of 
a 0  and  f;  obviously  failure  occurs  when 

t*  55  faQ  (B3) 

and  at  this  point  the  value  of  the  strain  is* 


Yf 


Yt 


2  1/3 

3(2-vs)  *  f ,  q0  . 

(l-vs)l/3  12GS2 


(B4) 


Substituting  the  properties  of  quartz  (Table  3)  in  eqn.  (B4)  we  obtain 
Yt  -  4.5  x  10-3  (0o)2/3  (85) 

with  oQ  in  psi. 


*  The  value  of  vs  ■  0.15  used  here  was  obtained  from  White  (1964)  and 
other  sources,  and  is  different  from  vs  =  0.31  used  in  a  previous 
report  (Dobry  et  al.,  1982).  The  value  vs  =  0.15  is  more 
representative  of  quartz;  as  a  result,  the  values  of  the  threshold 
strain  Yf  computed  using  Eq.  A66  and  vs  =  0.15  are  slightly  different 
from  those  originally  obtained  by  Dobry  et  al.  (1982). 
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B.2  The  Body  Centered  Cubic  Array  Subjected  to  Trlaxlal  Loading 
Consider  the  bcc  array  shown  in  Fig.  25,  initially  subjected  to  an 
isotropic  stress  o0  and  consider  an  additional  force  Pa  acting  in  the  X2 
direction.  Fig.  Al).  The  bcc  array  is  subjected  to  an  initial  isotropic 
stress  o0  and  the  stress  <*22  is  increased  monotonlcally  from  o0  to  a 


value  o0+oa,  at 

which  sliding  occurs  in  the  array. 

In  this  case,  the 

loading  path  is 

as  follows 

at  t  *  tQ 

P11  =  p22  =  p33  =  po 

(B6) 

at  t  ■  tj 

pll  "  p33  =  po 

(B7) 

p22  "  po  +  pa 

(B8 ) 

The  contact  forces  are  (eqns,  A5,  A6,  A7) 


VT~ 


dN3!  *  *  ~  ^PU  "  dl?33^  *  ° 


(B9) 


dN32  =*  ±  ^4  [dPU  "  2d?2  +  dP3 3 ]  *  ~[2  dPa 


(BIO) 


dN33  *  ±  [dPll  +  dP22  +  dP33]  =  ~J2  dPa 


(Bll) 


and  the  ratio  between  the  increment  of  the  tangential  force  and  the 
increment  of  the  normal  force,  8»  is 


«  ,  ,  ^21  ,  vO" 

6  dN  d^33 


(B 1 2) 


In  this  case  B  >  f»  the  coefficient  of  interparticle  friction,  therefore 
the  values  of  the  compliances  are:  (Mindlin  and  Deresiewicz,  1953) 
a)  Normal  Compliance,  Cn 

i  V  Q  /n  1  \ 


2G, 


’sa 


»  *  j  m 


(B 1 3) 


this  expression  for  being  valid  no  matter  what  the  loading  history  of 
the  spheres  is.  Now 


a  =  aQ  (1  +  0L) 


1/3 


(B14) 


where 


3(l-vs)  3(l-vs)  3 

8  R  =  -  —  P„  *  R 


8Gs 


8Ge 


(B13) 


Also  6  »  — 


T 

f 


The  vertical  displacement,  6^,  has  two  components;  a  6^  and  a  6^ 


component: 


2  *  2  • 
^22  =  533  +  —  623 

J  O 


(Bi6) 


the  6^^  component  is  computed  as  follows: 


8633  l-vs  ,  .-1/3 

C„  -  — ~  =  — — —  (1  +  0L)  u  J 


"n  dN  2G, 


(B17) 
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f inally 


( B20 ) 


633  =  2VT  R(1-vs)2/3  (■ 


1/3 


4V3Gq2 
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[(1 


.  0  2/3 

)  -  1] 

T  rt  '  J 


3  °c 


The  6^  component  will  be  computed  from  the  tangential  compliance  since 
it  is  a  tangential  displacement;  the  tangential  compliance  is  (Mindlin 
and  Deresiewicz,  1953): 


2~vs 

4Gsa 


[9  +  (1  -  0 ) f  1  - 
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-1/3 
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(B21) 


Simplifying: 
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now 
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ct  e[i  +  ef) 


-1/3 


+  ct  ( 1-9 )[  1  +  ~  (9-l)rI/3 

o  Lo 


(B24) 


where 
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Now 
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finally 
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[f1  +  (7  -0  -  1]}  <B26> 

8  r  o 


Simplifying  further  and  expressing  the  displacement  in  terms  of  stresses 
we  obtain: 


62 3-V3f R 


2-Vc 
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(l-vs)l/3  4VTG  2 
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2/3 


u  38  o0J  J  L'  J3f  ao 


2/3 


-1 


(B27) 

In  order  to  compute  the  vertical  displacement  of  the  array,  622 *  we  have 
to  substitute  equations  (B20)  and  (B27)  into  eqn.  ( B 1 6 ) : 


2  .  2 

622  =  ~  633  +  —  623 

VT  Ve 


(B28) 
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(B29) 


From  the  above  equation,  (B29),  the  stress-strain  relation  may  be  com¬ 
puted  for  the  bcc  array  under  isotropic  loading,  and  it  is: 


:22«( - — )  {V5Il-vs)"J  (I+-— )  -1  +—  f  - — 

4VTGs2J  3»oJ  J  4  ( l-vs )  1/3 


*  [f(l+^"—)  -  1]  [(1  +  (~  -  D^—)  -  1]]}  ( B 30 ) 

L  38  ao  8  3f  ao 


The  above  equation  is  plotted  in  Fig.  3/a  for  various  values  of  a0  and  f. 


At  this  point,  the  value  of  aa/a0  which  causes  failure  in  the  array  must 


be  determined.  Failure  is  defined  as  sliding  at  the  contacts;  this  time, 


since  the  bcc  array  is  statically  determinate,  failure  in  one  contact 


implies  failure  of  the  array.  Furthermore,  because  of  the  symmetry  of 


the  array,  failure  will  occur  simultaneously  at  all  contacts.  Sliding 


will  occur  when 


(B31) 


We  know  that 


T  ^A 

T  =  -  A  cr¬ 


am) 


M  ^3~" . 

N  "  —  A  a- 


(B33) 


VT" 

NQ  =  —  A  ar 
4  c 


(B34) 


where  A  is  the  area  of  the  face  of  the  bcc  array.  Then: 
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T  -  f(N0  +  N) 

VT\  ,VT  .  .  >/3  . 

—  A  o00  *  f( —  A  a-  +  —  A  o0) 
12  22  12  4  0 


“  °22  "  f<4  °  +  °  5 

3  22  3  a  o 


122  ,/?  .  f  i  ,  f 

°o  L  3  3 J 


and  finally,  aa/oQ  at  failure,  (oa/a0)  is 


f  2-f 


If  022  **  °o  +  aa  (total  stress),  then 


r°22.  f  ,  . 

'  zo 

f  3  3 


(B35) 


(B36) 


(B37) 


(B38) 


(B39) 


(B40) 


in  terms  of  total  stress  (Fig.  24) 

In  order  to  compute  the  strain  at  failure,  e22^»  we  most  substitute 
the  equation  for  (022/^0^ »  e9n*  ( B 3 9 )  into  the  stress-strain  relation, 
eqn.  (B30).  Doing  this  we  obtain 


,  °°  >  ,  .2/3  rn  ,  f  ,  n  .  3^2  2-vs 


s  * 


f  */  J  c  pr  1 

[1(1  +— — )  -u-u  +  (-  -  1)  -  1)]} 


(B41) 
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for  the  properties  of  quartz  (Lambe  and  Whitman,  1964,  Ko  and  Scott,  1967, 
White,  1964) 

Gs  =  4.783  x  106  psi 

vs  =  0.15 

f  =  0.5 

the  strain  at  failure  is 

e22^  =  3.438  x  10-3  a0^/3  (in  percent)  (B42) 


and  Yt  =  e22f  for  vs  =  °» 

which  is  the  threshold  strain  yt  for  the  array.  This  equation  is  plotted 
in  Fig.  27  together  with  the  other  expressions  for  yt  for  the  other 
arrays  for  easy  comparison. 

B.3  The  Face  Centered  Cubic  Array  Subjected  to  Triaxial  Loading 
The  triaxial  loading  of  a  Face  Centered  Cubic  Array  was  solved  by 
Brauns  and  Leussink  (1970).  In  this  work,  the  stress-strain  relationship 
is  not  obtained  for  the  whole  range  of  values  of  022/^0  But  only  for 
those  which  make  the  array  statically  determinate.  It  is  extremely  hard 
to  determine  the  values  of  compliances  for  cross  anisotropic  loading 
(Duffy  and  Mindlin,  1957);  therefore  once  sliding  occurs  and  the  number 
of  contacts  decreases  from  12  to  8,  the  array  becomes  statically  deter¬ 
minate  and  it  is  possible  to  compute  a  stress-strain  relation  up  to  the 
point  which  the  array  fails  (range  b  -Fig.  B2).  Once  the  array  has 
failed,  simple  geometric  considerations  make  the  computation  over  the 


strain  range  c  possible  (Fig.  B2). 

Consider  the  fee  array  in  Fig.  Bla  and  the  cross  anisotropic  loading 
as  shown.  The  free  body  diagram  of  the  octant  of  sphere  A  is  shown  in 
Fig.  Bib.  The  equilibrium  conditions  yield 


N  +  T  =  yfl  R2  oi 


N  -  T  +  NX  =  2>/2~  R2  03 


(B43) 


(B44) 


Also,  the  sum  of  the  displacements  around  a  closed  path  must  vanish  (Duffy 


and  Mindlin,  1957),  which  yields 


a^-a+6-0 


The  normal  compliance  is  found  to  be 


(B45) 


[3(1-vs)Gs2RN] 1/3 


(B46) 


and  the  Tangential  Compliance 


ct 

dT 


The  strain  is 


2[3(1-vs)Cs2RN]1/3 


*  [f-» 

L  JT 


i-f5 

dT 


t  1/3 
[1  -  —1 
fN 


(B47) 


ei  1  =  —  (a  +  6) 
11  2R 


(B48) 
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Integrating  the  compliances  and  substituting  the  results  into  eqn.  (B48), 
we  find  after  transforming  the  forces  into  stresses  that  the  strain,  ejj, 
is  given  by 


„  r3/2(l-us)2/3  l  ,qj_  1 27 3  1-f  ,,  2/3 

11  L  §g^  T  1  L  i+f  ^033^  1  1+f  ^033'^  33 


(B49) 

this  expression  being  valid  only  for  the  range  of  (013/033)  in  which  the 
array  is  statically  determinate,  that  is  8  contacts  per  sphere  (range  b 
in  Fig.  Bl). 

Failure  is  defined  again  as  sliding,  but  this  time  when  the  whole 
array  fails,  that  is  when  the  number  of  contacts  from  8  reduces  to  none. 

Using  the  same  criteria  as  in  the  other  arrays,  the  critical  stress 
ratio  at  which  the  array  fails  is 


(B30) 


The  above  equation  is  plotted  in  Fig.  31b 


Figure  Bl.  (a)  Representative  Unit  Volume  of  an  fee  Array  with 
State  of  Stress  (b)  Sphere  Centered  a  Apex  A  with 
Applied  Stresses,  Contact  Forces  and  Displacements 
(Brauns  and  Leussink,  1970)  . 


Figure  B2.  Triaxial  Compression  of  an  fee  Array  of  Glass 
Spheres:  Analytical  and  Experimental  Results 
(Brauns  and  Leussink,  1970)  . 
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C  SUBROUTINE  UMAT  TO  88  USED  WITH  TH8  FIMITI  ELEMENT  PROGRAM 
C  ABAQOS 

C  . 

C 

S0BROOTI HE  OMAT ( I , STRESS , DDSDOE , STRAN , DSTRAM , NDI , NSHR , NTENS , M ) 
IMPLICIT  REAL*8(A-H,0-Z) 

C0M40N/BMATRL/GI ,  2 ,  ON ,  SR ,  DEBOG 
COMMON/BHI ST/ AHI ST ( 1,3004) 

LOGICAL  FLAG 1,PLAG2, DEBOG 

DIMENSION  STRESS (NTENS) , DDSDOE (NTENS, NTENS ), STRAN (NTENS) 
DIMENSION  DSTRAN ( NTENS ) ,DSIGMA(4) 

AREA-4. *SR*SR 
DO  141  JI-1,4 
DO  142  J2-1 , 4 
DDSDOE ( J 1 ,J2)-0.D0 

1 42  CONTINUE 
141  CONTINUE 

C 

C  CHANGE  COMPRESSION  INTO  POSITIVE 
C 

DO  143  J3-1 ,3 
STRESS(J3 ) --STRESS (J3 ) 

STRAN ( J 3 ) --STRAN ( J 3 ) 

DSTRAN  ( J  3 )  —DSTRAN  ( J  3  > 

143  CONTINUE 

C 

C  TRANSFORM  STRESSES  TO  FORCES 
C 

TN 1 1 -STRESS ( 1 ) * AREA 
7X2 2 •STRESS ( 2 ) ‘AREA 
TN3 3-STRESS ( 3 ) ‘AREA 
T1  2-STRESS ( 4 ) "AREA 
C 

C  COMPUTE  CURRENT  TOTAL  NORMAL  DISPLACEMENTS  FROM  STRAINS 
C 

DEL 1 1 • ( STRAN ( 1 ) ♦DSTRAN ( 1 ) )*2.*SR 
DEL22- ( STRAN ( 2 ) +DSTRAN ( 2 ) ) *2 . *SR 
C 

C  COMPUTE  THE  INCREMENTAL  SHEAR  DISPLACEMENT  FROM  STRAINS 
C 

DOEL 1 2- 2 . ‘DSTRAN ( 4 ) *SR 
C 

C  COMPUTE  NORMAL  FORCES  AND  INCREMENTS 
C 

CALL  NORMAD(DEL1 1 , TNI  1 ,TNNA,DTNA) 

CALL  NORMAD(DEL22,TN22,TNNB,DTNB) 

C 

C  TRANSFORM  FORCES  TO  STRESSES 
C 

OS I GMA ( 1 ) -DTNA/AREA 
DSIGMA( 2 ) -DTNB/AREA 
DSIGMA(3)-0. 

C 

TOL-1 .E-16 
EPSL-t .E-16 
ITMAX-500 
DETA-0.D0 


c  & 


*,  »l,  I'-h*.  IVi'tit.  »•.  !*.  |*.  It. 


DETB-O.DO 

DTYA-0.00 

101-1 

102-1503 

C 

C  COMPOTE  SEED  FOR  THE  ITERATION, 

C 

IP  (DSTRANU)  .BQ.  0.)  THEM 
DEXA-0. 

DEXB-0. 

ITER-0 

DTXA-0. 

oTsa-o. 

GO  TO  1167 
ENDIP 
C 

IF ( (TNNA-TNNB)  .EQ.  0.)  THEN 
WRI TE ( 3 , 1 37 5 ) TNNA , TNNB 
1375  FORMAT ('TNNA-TNNB* ,1*,P10.6,2X,P10.6) 

DEXA-DDBL 12/2.00 
DEXB-DDEL 12/2.00 
GO  TO  1169 
ELSE 

□EZTA-DDEL 12/2.00 
DEXTB-DOEL 12/2.00 
ENDIP 

IP ((TNNA-TNNB)  .GT.  0.)  GOTO  1311 
C 

C  BEGIN  ITERATION  B 
C 

CALL  CONPLD ( M, I 02 , TNNB , DTNS , DEXTB , DETB , OTXB , DTTT , I DEXT ) 

C 

C  GIVEN  THE  SHEAR  STRAIN  ITERATE  TO  PINO  THE  COMMON  SHEAR  FORCE 
C 

DO  1 166  ITER-1 ,ITMAX 
DTTB-O.DO 

CALL  CONPLF ( M , I D2 , TNNB , DTNB , DTXB , DTYB , DEXB , DEYB , I DEXB ) 
DEXA-DDEL1 2-OEXB 

CALL  CONPLD(M, ID1 , TNNA , DTNA , DEXA , DEYA , DTXA , DTYA , I DEXA ) 

I F ( DABS ( DTXA-DTXB )  .LT.  TOL)  GO  TO  1167 
DTXB- ( DTXA+DTXB ) /2 . DO 
1166  CONTINUE 
C 

IF (  ITER  .GE.  ITMAX  .AND.  IDEXA  .EQ.  1)  THEN 
DTXB-0.0 
DTXA-0. 0 
OEXB-O.O 
DEXA-DDEL1 2 
I  TER- ITMAX 
GO  TO  1167 
ENDIP 
C 

IF  (ITER  .GE.  ITMAX  .AND.  I DEXB  .EQ.  1)  THEN 
DTXB-0.0 
DTXA-0. 0 
DEXB-DDEL1 2 
DEXA-0. 0 
ITER-ITMAX 
GO  TO  11 67 
ENDIP 


'SWV* 


wv  "■r  ■-■  <>  v  \>*  "jw^  f9^y^s^ww^yf 

>M 
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WRITE(3,1954)I 

1954  FORMAT ( ' WELL  ITER  WAS  I THAI  ROT  IDE!  WAS  NOT’, 1!, 15) 

C 

1311  CONTI HUE 
C 

C  BEGIN  ITERATION  A 
C 

CALL  C0NPLD(M,ID1 , TNNA , DTNA , DEXTA , DETA , DTXA , DTT , I DET ) 

C 

C  GIVEN  THE  SHEAR  STRAIN  ITERATE  TO  FIND  THE  C0M40N  SHEAR  FORCE 
C 

DO  1168  ITBR-1 ,ITMAX 
DTT  A  >0.  DO 

CALL  CONPLF(M,ID1 , TNNA , DTNA , DTXA , DTTA , DSZA , DETA , I DEXA ) 
DEXB-DDEL12-D8XA 

CALL  CONPLD ( M , I D2 , TNNB , DTNB , DEXB , DETB , DTXB , DTTB , I DEXB ) 

I F ( DABS ( DTXA-DTXB )  .LT.  TOL)  GO  TO  1167 
DTXA-  ( DTXA  DTXB )  /2 .  DO 

1168  CONTINUE 
C 

I F (  ITER  .GE.  ITMAX  .AND.  IDEXA  .EQ.  1)  THEN 
DTXB-0.0 
DTXA-0.0 
DEXB- 0.0 
DEXA-DDEL 1 2 
I TER- ITMAX 
GO  TO  1167 
ENDIF 
C 

IF  (ITER  .GE.  ITMAX  .AND.  I DEXB  .BQ.  1)  THEN 
DTXB-0.0 
DTXA-0.0 
DEXB-DDEL1 2 
DEXA- 0.0 
ITER-ITHAX 
GO  TO  1167 
ENDIF 
C 

1 167  CONTINUE 

WRITE( 1,111 1 )I ,ITER, DEXA. DEXB, DTXA, DTXB 
1111  FORMAT (215, 1X,F10.8,1X,F10.8,1X,F10.8,1X,F10.8) 

C 

C  UPDATE  CONTACT  HISTORY  AND  STORE  INTO  AHIST 
C 

1169  CONTINUE 

CALL  CONPD(M,ID1 , TNNA , DTNA , DEXA , DEYA , DTXA , DTYA , I DEX 1 ) 

CALL  CONPD ( M , I D2 , TNNB , DTNB , DEXB , DETB , DTXB , DTYB , I DEX2 ) 

C 

DTXM- ( DTXA*DTXB ) /2 . DO 
DS I GMA ( 4 ) • DTXM/ AREA 
WRITE(6, 1 863 ) DSIGMA ( 4 ) 

1863  FORMAT(F20. 10) 

C 

C  COMPUTATION  OF  THE  INCREMENTAL  STIFFNESS  MATRIX,  STIFF 
C 

AB-3 . DO* ( 1 . DO-UN) /( 8 .D0*GI ) 

ABA-(AB*TNNA*SR)**( 1 ,/3. ) 

ABB- ( AB*TNNB*SR) **( 1 ./3. ) 

DDSDDE( 1 , 1 ) -2 . D0*GI * ( AB*TNNA*SR ) ** ( 1 ./3. )/( 1 .DO-UN) 
DDSDDE(2,2)-2.D0*GI*(AB*TNNB*SR)**( 1 ./3. >/( 1 .DO-UN) 


onnvo  n  n  o  o  non  onn  no  o  non  non 


DDSDDB ( 3 , 3 )  *  1 .DO 
IP  (DSTRANU)  .BQ.  0.)  THEN 
DDSDDB(4,4)-1 .DO 
ELSE 

DDSDDB  (4,4)  -DS I GMA (  4 )  /DSTRAN  ( 4 ) 

END  IP 

IF  ( DDSDDE (4.4)  .LB.  EPSL)  DDSDDB ( 4,4) -EPSL 

UPDATE  STRESS  VECTOR 

STRESS  ( 1 ) -STRESS  ( 1  )+DSIGMA<  1 ) 

STRESS ( 2 ) -STRESS ( 2 ) +DSIGMA ( 2 ) 

STRESS(3)-0.D0 

STRESS  ( 4 )  -STRESS  ( 4 )  +DS I GMA  ( 4 ) 

CHANGE  COMPRESSION  INTO  NEGATIVE 

STRESS ( 1 ) -“STRESS ( 1 ) 

STRESS ( 2 ) --STRESS ( 2 ) 

RETURN 

END 


SUBROUTINE  CONPD(M,IAD1 , TNN , DTN , DEE , DET , DTX , DTT , I DEZ ) 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/ BMATRL/G I , Z , UN , SR , DEBUG 
COMMON/BHI ST/ AHI ST ( 1 , 3004) 

DIMENSION  XKKI  (100),XKKS  (100), HHI  (100), HHS  (100), WZ  (100), 
$WT( 100), TTXS< 100) .TTYS  (100),TTXI  (100),TTYI  ( 1 00 ) , AAXS( 1 00 ) , 
SAAYS  (100),AAXI  (100),AAYI  (100), IC  (100) 

LOGICAL  FLAG 1.PLAG2, DEBUG 

- INITIALIZATION - 


I DEI-0 

IN- I ADI 

DEBUG-. FALSE. 

NG"IDINT(AHIST(M, IAD1 ) ) 

DTX-0. 

DTT-0 . 

NGROUP- 1 00 

PREVIOUS  NORMAL  FORCE,  TNI 


TNI -TNN-DTN 

MAKE  SURE  THERE  IS  A  CONTACT;  FOR  THIS,  THE  NORMAL  FORCE  HAS  TO 
BE  POSITIVE--IDEX-2  MEANS  THAT  THERE  IS  NO  CONTACT 

IF  (TNN.GT.0)  GO  TO  9827 

IDEX-2 

RETURN 

827  IF(NG.NE.O)  GOTO  9820 

CALL  THE  ELASTIC  MODULUS  OF  THE  CURRENT  NORMAL  FORCE 

CALL  EMOD(TNN.HO) 

TX«H0*DEX 
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TM-DSQRT ( TX*TX+TT*TT ) 
OTH-TM 

IP(DTM.BQ.O. )  GOTO  9840 
IP  (TM  .GE.  Z*TNN)  THEM 

IDEX-1 

TX-TX*Z*TNN/TM 

TY-TY*Z*TNN/TM 

TM-Z*THM 

END  IP 

DTX-TX 

DTT-TT 


C 

C -  INITIALIZE  THE  NEW  CONTACT  HISTORY 

C 

9840  NG-1 
C 

C  INITIAL  GROUP  CHARACTERISTICS 

C  - MODULI  AND  SIZES - 

C 


HHS ( 1 ) »H0 
HHI(1 )-0. 

XKKS ( 1 )-Z*TNN 
XKKI ( 1 )«0. 


C 

C  -  INITIAL  SURFACE  POSITIONS  - 

C  COORDINATES  OP  CENTRES  OP  MAX  BOUNDARY  YIELD  CIRCLE 

C 


AAXS ( 1 ) "0 . 

AAYS(1)«0. 

C 

C  COORDI ANATES  OP  CENTRES  OP  MIN  BOUNDARY  YIELD  CIRCLE 

C 

AAXKO-TX 
AAYI ( 1 ) -TY 
C 

C - INITIALS  REFERENCE  POINTS - 

C 

TTXS<1)-0. 

TTYS( 1 )-0. 

TTXI ( 1 ) »TX 
TTYI ( t )-TY 


C 

C - KIND  OP  INITIAL  GROUP 

C 

IF(IDEX  .EQ.  0)  THEN 
C 

C - NO  SLIDE - 

C 


IC( 1 )-0 
WX(  1  )  - 1  . 

WY (  1  )-1  . 

VX-1  . 

VY»1  . 

GO  TO  9910 
ELSE 
C 

C -  SLIDE  OCCURS 

C 

I C  (  1  )  -  1 
WX(  1  )-TX/TM 


AiV.V.V.V.V.V  Vi 
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WT(1  J-TT/TM 
TTXS ( 1 )-TX 
TTYS(1)-TY 
IDBX-1 
GO  TO  9910 
END  IF 
C 

C -  BEGIN  CALCULATION  WITH  THE  OLD  CONTACT 

C 

9820  DEM-OSQRT(DEX*DEX+DEY*DEY) 

C 

C —  CHECK  IF  THE  MEMORY  IS  ENOUGH 
C 

IF(NG.LT.NGROUP)  GO  TO  9790 

IOEX-4 

RETURN 

C 

C  -  GET  THE  HISTORY  OF  THE  CONTACT  FROM  THE  MEMORY 
C 

9790  CONTINUE 

DDTN-AHIST(M,IAD1-M  ) 

DO  9810  I ■ 1 , NG 

IAD2"IAD1+I+1 

XKKI (I ) -AH 1ST (M, I AD 2 ) 

XKKS ( I ) - AHI ST ( M , I AD2+NGROUP ) 

HHI (I )-AHIST(M,IAD2*2*NGROUP) 

HHS(I )«AHIST(M,IAD2*3*NGROUP) 

Will  )-AHIST(M,IAD2«-4*NGROUP) 

WY ( I ) -AHI ST ( M ,  I AD2* 5*NGROUP ) 

TTXS ( I ) -AHI ST ( M , I A02+6*NGR0UP ) 

TTYS ( I ) - AHI ST ( M , I AD2+ 7 *NGROUP ) 

TTXI ( I ) -AHI ST ( M , I AD2 >8 *NGROUP ) 
TTYI(I)-AHIST(M,IAD2*9*NGROUP) 

AAXS ( I ) «AHIST(M, IAD2+ 1 0*NGROUP) 
AAYS(I)-AHIST(M,IAD2-M  1*NGROUP) 

AAXI ( I ) - AH I  ST ( M , I AD 2 *1 2 *NGROUP ) 

AAYI ( I ) -AHIST(M, IAD2* 1 3*NGROUP) 

9810  IC(I)-IDINT(AHIST(M,IAD2-M4*NGROUP) ) 

TX-TTXI (NG) 

TY-TTYI (NG) 

C 


C  - END  OF  INITIAL  REFERENCE  POINTS- 

C 

IF(DEM.EQ.0.)  GOTO  9900 
C 

C - UPDATING  PARAMETER  BEFORE  INCREMENT 

C 


HSL-HHS(NG) 

HIL-HHI (NG) 

IF ( HIL . LE . 0 . )  HIL-HSL 
XKL-XKKS(NG) 

IF(TN1 .GT.0. )  CALL  ELAD ( XKL , TN 1 , EML ) 
HML-0.D0 

IF(EML.NE.0. )  HML-XKL/EML 
NGL-NG 
C 

C -  CALCULATE  THE  CORRECT  FORCE  INCREMENT 

C 

DTX»HSL*DEX 

DTY-HSL*DEY 


no  non  n  no  non  ua  non  onono 


v*vvvnrvTrjWwwv,’^.'^"A.V‘  ~j>  ->  .-v vwv^vvs,  -.-  *.-  s-V-' 


OTH-HSL*DBM 


-UNIT  NORMAL  VECTOR- 


IFdC(NGL)  .EQ.O  .OR.  DTN  .GT.  0.)  THEN 

VX-DEX/DEM 

VY-DEY/DEM 

ELSE 

VX-WX(NGL) 

VY-WY(NGL) 

VXL-VX 
VYL-VT 
END  IP 

DBN-DTX*VX*DTY*VY 
DBT-DABS ( DTX*VY-DTY*VX ) 

FIND  THE  CORRESPONDING  PORCE 

I P  < DBN . LE . Z*DTN )  GOTO  9900 
IFdC(NGL)  .EQ.O)  THEN 
DTM-(DEM*2*DTN*( 1 ./HIL-1 ./HSL) ) *HIL 
ELSE  I P ( HML .EQ.O.)  THEN 

DTM-DEM*(DEM+Z*DTN*( 1 ./HIL-1 ./HSL) ) *HIL/(DEX*VX+DEY*VY) 
ELSE 

DTM-DEM*(DEM+Z*DTN*( 1 ./HIL- 1 ./HSL) )/( (DEX*VX+DEY*VY)/HIL 
$  ♦DABS (DEX*VY-DEY*VX ) /HML) 

ENDIP 

DTX-DEX*DABS ( DTM) /DEM 
DTY-DEY*DABS ( DTM ) /DEM 
DBN-DTX*VX*DTY*VY 
DBT-DABS ( DTX*VY-DTY*VX ) 

900  TN-TNN 

TX-TX*DTX 

TY-TY+DTY 

TM-DSQRT ( TX*TX*TY*TY ) 


-UPDATING  FORCE  INCREMENT- 


TNL-TN 

TXL-TX 

TYL-TY 


-CURRENT  ELASTIC  MODULUS- 


CALL  EMOD(TN.HO) 


-CHECK  IF  THERE  IS  FAILURE- 


IF(TM.LT.Z*TN)  GO  TO  800 
IDEX-1 

TX»TX*Z*TN/TM 

TY-TY*Z*TN/TM 

TM»Z*TN 

DTX-TX-TTXI (NG) 
DTY-TY-TTYI (NG) 

GOTO  9840 


-UPDATE  THE  SIZES- 


m 


& 
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800  00  211  IG-1 ,  NG 


C 

C- 

C 


RADIOS  OF  YIELD  CIRCLES  OF  GROUP  (IG) 


XXXI (IG) -XXXI (IG)+Z*DTN 
XXX S (IG) -XXXS( IG) *Z*DTN 


C-—  THIS  IS  TO  MAXS  SORB  THAT  HO  MODOLUS  IS  LARGER  THAH  HO  — 
C 

IF(XKXS(IG) . LT.O. }  XKKS(IG)-0. 

211  IF<XKXI (IG).LT.O.)  XXXI(IG)«0. 

C 

C - FIX  THE  SIZE  (RADIOS)  OF  THE  FAILURE  SORFACE - 

C 

XKKS ( 1 ) -Z*TN 


C 

C - UPDATING  THE  MODULI - 

DO  215  IG- 1 ,NG 

HHI (IG)-( ( 1 .-(XXRS(IG))/(Z*TN))**( 1 ./3. ) )*H0 
HHS(IG)-( ( 1 . -(XKKI (IG) )/(Z*TN) )**( 1 ,/3. ) )*H0 
215  CONTINUE 

C  G-- (TTYI ( NGL ) -TTYS ( NGL ) ) *DTX* ( TTXI ( NGL ) -TTXS ( NGL ) ) •DTY 

G--DDTN*DTX+ ( TTXI ( NGL ) -TTXS ( NGL ) ) *DTN 


C 

C - UPDATE  REFERENCE  POINT  WHEN  IC-1 - 

C 

DO  213  IG-1 ,NG 

IF(IC(IG) . EQ  0)  GO  TO  213 

TTXS  (I  G ) -TTXS  (I G )  f  Z*DTN*  WX  (I G ) 

TTYS  ( I G )  -TTYS  (I  G )  ♦  Z *DTN*WY  (I G ) 

TTXI  (IG) -TTXI  <IG)fZ*DTN*WX(IG) 

TTYI  (IG ) -TTYI  (IG ) *Z*DTN«WY (IG ) 

213  CONTINUE 

IF  (DTN)  849,  850,  851 

849  IF  (TX  .EQ.  0.  .AND.  TY  .EQ.  0.  .AND.  DTK  .EQ.  0. 
$ .AND.  IC(NG)  .EQ,  0.  .AND.  XXXS(NG)  .GT.  0.)  THEN 

XK-XKKS(NG) 

H-H0* ( 1 -XX/ ( Z*TN ))**(l./3.) 

AX-AAXS ( NG ) 

AY-AAYS(NG) 

ICI-0 
GOTO  60 
END  IF 
GO  TO  272 

850  IF  (DEM  .EQ.  0 . ) RETURN 
GO  TO  272 

851  IF  (ABS(DBN)  .LT.  Z*DTN)  THEN 

IF  (IC(NG)  .EQ.  0.  .AND.  G  .EQ.  0.)  THEN 
XK-XKKS ( NG ) 

H-H0M 1 .-XK/(Z*TN) )**( 1 ,/3 . ) 

AX-AAXS (NG) 

AY-AAYS(NG) 

ICI-0 
GO  TO  60 
END  IF 
GO  TO  273 
END  IF 

IF  ( IC(NG)  .EQ.  1  .AND.  DBN  .EQ.  Z*DTN)  THEN 
XK-XKKS (NG) 


t»  » 

>:v?! 


W 


non  n  n  n  o  o  n  n  -  k>  ro  o  o  n  n  non 


H-H0*(1.-XR/(Z*TM))**(1./3.) 

AX-AAXS(NG) 

AY-AAYS(NG) 

ICI-1 
GO  TO  60 

BMOir 

GO  TO  272 

- A  HEW  ELASTIC  GROUP  IS  CREATED' 

- "NG’IS  INCREASED  TO  "NG*! . 


273  CONTINUE 
NG-NG-H 
XK-X*DTN 

H- ( ( t . -DTH/TN ) ** ( 1 . /3 . ) ) *H0 
AZ-AAZKNG-1  ) 

AY-AAYKNG-1 ) 

ICI-0 
GO  TO  60 

- TO  WHICH  REFERENCE  GROUP  DOES  THE  rORCE  POINT  BELONG - 

-  THIS  IS  THE  SUBROUTINE  SEARCH  OF  THE  ORIGINAL  PROGRAM  — 

72  IG-0 
NG-NGL 
76  IG-IG+1 

IF(IG.GT.NG)GOTO  273 
TAXS-TX-AAXS(IG) 

TAYS-TY-AAYS(IG) 

RS-DSQRT(TAXS*TAXS*TAYS*TAYS) 

IF(RS.GT.XKKS(IG) )  GO  TO  99 
TAII-TX-AAXI (IG) 

TAYI-TY-AAYI(IG) 

RI “DSQRT ( TAXI *TAXI ♦TAYI *TAYI ) 

IF (DEBUG)  WRITE{ 11,1234)  IG,TX,TY, RS ,XKKS ( IG) , RI , XKKI ( IG ) 
234  FORMAT( 2X, ' IG' , 51, ' TX' , lOX.’TT' , I  OX, 'RS' , 9X, ' XXKS ' , 9Xf 
. 'RI' ,5X, ’XKKI’ ./I5.6E12.4) 

IF(XRRKIG).GT.RI)  GO  TO  276 

- APPROXIMATE  LOADING  MODULUS  AND  NUMBER  OF - 

- YIELD  SURFACES  AFTER  A  LOAD  INCREMENT - 


99  CONTINUE 
NG*IG 

I F ( DTN . LT . 0 . AND . NG . NE . 1 . AND. IC(NGL) . EQ. 0 . AND.G .EQ. 0 . )GOTO  9910 
I F ( I C ( NG ) . EQ . 1 )  THEN 

-  raE  FORCE  POINT  FALLS  INSIDE  A  GROUP  OF  THE  SECOND  KIND 


XK-0. 

TEMPI  1--2.*(  ( TX-TTXS  ( NG )  )  *WX  ( NG )  ♦  ( TY-TTYS  ( NG )  )*WY(NG)  ) 
I F ( TEMP 1 1 . NE. 0 . )  XK-( (TX-TTXS ( NG ))*( TX-TTXS (NG ) ) 

$  ♦ (TY-TTYS (NG) )* (TY-TTYS (NG) ) )/TEMP1 1 

ELSE 


-— -  THE  FORCE  POINT  FALLS  INSIDE  A  GROUP  OF  THE  FIRST  KIND 

AA 1 ■ (TX-TTXS (NG) ) * (TX-TTXS (NG ) ) ♦ (TY-TTYS (NG ) ) * ( TY-TTYS ( NG )  ) 
AA2- ( TTXS ( NG ) -TTXI ( NG ) ) * ( TX-TTXS ( NG ) ) 

S  ♦ ( TTYS ( NG ) -TTYI ( NG ))*( TY-TTYS ( NG) ) 

AA3 ■ ( TTXS ( NG ) -TTXI (NG ))* (TTXS ( NG ) -TTXI ( NG ) ) 


n  n  n  nno  non  nnn  non 


$  +(TTYS (NG)-TTYI (NG) )* (TTTS (NG)-TTTI (NG) ) 

ALPH-AA3/ ( ZKRS ( NG ) *XKKS ( MO ) ) - 1 • 

BETA- ( AA3+AA2 ) /XKKS ( NG) 

GAMA-AA1 +AA2+AA2+AA3 
IF(ALPH.NE.O. )  GOTO  900 


ZK-0. 

IP(BETA.NE.O. )  XK -GAMA/ (2.* BETA) 

GOTO  901 

900  XK- <  BBTA-DSQRT ( BETA* BETA- ALPH*GAMA ) ) /ALPH 
EMDIf 

901  If ( .NOT. DEBOG)  GOTO  902 

WRITE( 11,134)  XK , TTXS ( NG ) , TTTS ( NG ) ,TTXI (NG) ,TTYI (MG) , 

$  WX(NG) ,  WT(NG) 

134  FORMAT(5X, 'XK’ ,9X, 'TTXS' ,8X, 'TTYS' , 8X , ' TTXI * ,  8X , ' TTTI ' , 

$  8X, '  WX*  ,8X,  ’WT'  ./7E12.5) 

If (IC(NG) .EQ.0. )  WHITE( 11,1 35)  AA1 ,AA2,AA3,ALPH, BETA, GAMA 

135  FORMAT( 5X, ' AA1 ' , 9X, ' AA2 ' , 9X, ' AA3 ' ,8X,'ALPH' , 8X , ' BETA ' , 8X , ' GAMA ' 

$  /6E12.5) 

902  IP(XK.GE.O .ANO.XX.LE. Z*TN)  GOTO  119 

WRITE ( 11,134)  XK , TTXS ( NG ) , TTYS(NG) ,TTXI (NG) ,TTTI (MG) , 

$  WX(NG)  ,WT(NG) 

If ( IC(NG) .EQ. 0. )  WRITE( 1 1 , 135)  AA1 ,AA2, AA3,ALPH, BETA, GAMA 

IOEX-3 

RETURN 


FINAL  PARAMETERS 


119  ICI-1 

H-H0«( 1-XK/(Z*TN) )**( 1 ./3. ) 
If(IC(NG).EQ.1)  THEN 
AX- TTXS  ( NG )  -XK*WX  ( NG ) 

AY-TTYS  ( NG )  -XK*  WY  ( NG ) 

ELSE 

AA4- ( H*H*H-HHI ( NG ) *HHI ( NG ) *HHI ( NG ) ) 

S  /(H0*H0*H0-HHI (NG) *HHI (NG) *HHI (NG) ) 
AX-TTXS ( NG ) ♦ ( TTXI ( NG ) -TTXS ( NG ) ) * AA4 
AY-TTYS ( NG )♦( TTYI ( NG ) -TTYS ( NG )) *AA4 
END  IF 

FINAL  PARAMETERS 


- OUTWARD  NORMAL  UNIT  VECTOR  WHEN  XK-0. - 

IF(XK.EQ.0.)  THEN 
VX-WX(NG) 

VY-WY(NG) 

DTM-DSQRT ( DTX*DTX*DTY*DTY ) 

If  (TTXI(NG)  .EQ.  TX  .AND.  TTYI(NG)  .EQ.  TY  .OR.  DTM  .EQ.  0) 
$  GOTO  555 
VX-DTX/DTM 
VY-DTY/DTM 

WHEN  THE  LOADING  IS  ALONG  THE  CONE  AND  NORMAL  FORCE  DECREASE 
ELSE 


OUTWARD  NORMAL  UNIT  VECTOR  WHEN  R(»XK)  IS  DIF  FROM  0. 
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VX-(TX-AX)/XX 

VY-(TY-AY)/XK 

BNOir 

5SS  IF(DEBUG)WRITS(n,  I97)IN,VX,VY,XK,H 

1 97  FORMAT  ( 3X, ' IN'  ,61,'VX'  ,  10X,'VT'  ,  1  OX, '  XX'  ,8X, 1 H'  ,/I  5 , 4E1  2 .4) 


C 

c 

C - GROUP  NUMBER  OF  THE  CURRENTLY  ACTIVATED 

c - YIELD  SURFACE  IS  *NG-M  * - 

C 

NG-NG-H 

C 

C - CONTINUITY  OF  SIZE - 

C 

XKKKNG-1 )  -XK 
HHS (NG“ 1 )-H 
C 


C  COORDINATES  OF  CENTER  OF  YIELD  CIRCLE  OF  MIN  RADI US (INFERIOR) 
C 

aaxKng-i  )-ax 
AAYKNG-1  )"AY 
C 

C - UPDATE  FOR  THE  REFERENCE  POINT - 

C 

IF(IC(NG-1).EQ.1)  GO  TO  60 
TTXI (NG-1 )-AX 
TTYKNG-1  )-AY 
AAZKNG-1  )-AX 
AAYI(NG-1 )-AY 


C 

C - UPDATING  THE  CONFIGURATION  OF  THE 

c - YIELD  SURFACES - - 


C 

60  CONTINUE 
XKKS(NG)-XK 
XKKI (HG) "0 . 

HHS  ( NG )  “HO 
HHI (NG) "H 
AAXS(NG)-AZ 
AAYS ( NG ) »AY 
AAXI ( NG ) «TX 
AAYI ( NG ) »TT 
TTXS ( NG ) «TX 
TTYS ( NG ) «TY 
TTXI (NG) »TX 
TTY I (NG)-TY 
WX(NG)-VX 
WY  ( NG )  *VY 
IC(NG)-ICI 

IF(ICI.EQ.O)  TTXS (NG ) "AX 
IF(ICI.EQ.O)  TTYS(NG)-A 
9910  DO  109  IG-I.NG 

IF(DEBUG)WRITE( 13, 124)IN,IG,HHS(IG) ,HHI(IG) ,XKKS(IG) , 

1 XKKI ( IG) , IC( IG) 

109  IF ( DEBUG )WRITE( 15, 124) IN, IG,AAXS( IG) , AAYS ( IG ) , AAXI ( IG) , 
1AAYI (IG) ,IC(IG) 

IF ( DEBUG )WRITE( 13,7) 

IF ( DEBUG )WRITE( 15,7) 

124  FORMAT ( IX, 215, 4E 12. 4, 15) 

115  FORMAT ( IX, 215, 2F 12. 4, 15) 


ic'  4i 


’  li'li'  .4  ■41'aH‘A'alfeWi  -«»  .>r  u>  a»>  a-a'lj  < 


7  FORMAT (IX,1 


UPDATE  THE  HISTORY  BEFORE  LEAVING  THE  SUBROUTINE 


AHISTtM, I ADI )-DFLOAT(NG) 

AHISTtM, IAD1 +1 )-DTN 
DO  9960  1-1, NG 
IA02-IA01+I+1 
AHISTtM, IAD2 ) -XKKI ( I ) 

AHISTtM, I AD2+NGROUP ) -XKKS ( I ) 

AHISTtM, IAD2+2*NGROUP)-HHI (I ) 

AHISTtM, IAD2+3*NGROUP)-HHS( I ) 

AHISTtM,  IAD2+4*NGROUP)-WX(  I ) 

AHISTtM, IAD2«-5*NGROUP)-WY(I ) 

AHISTtM, IAD2+6*NGROUP)-TTXS( I) 

AHISTtM, I AD2+7*NGROUP) -TTYS ( I ) 

AHISTtM, IAD2+8*NGROUP) -TTXI ( I ) 

AHISTtM, IAD2*9*NGROUP) -TTTI ( I ) 

AHISTtM, IAD2+10*NGROUP)-AAXS(I ) 
AHISTtM, IAD2+1 1 *NGROUP ) -AAYS ( I ) 
AHISTtM, IAD2+1 2*NGROUP)-AAXI ( I ) 
AHISTtM, IAD2+1 3*NGROUP ) -AAYI ( I ) 

9960  AHISTtM, I AD2+ 1 4 *NGROUP ) -OFLOAT ( ICt I ) ) 
RETURN 

C - OTHER  FORMATS - 

333  FORMAT (I5,7E12.4) 

END 


SUBROUTINE  EMOO t TNP , HO ) 


MATERIAL  DEPENDENT  SUBROUTINE  COMPUTES  THE  ELASTIC 
MODULUS  CORRESPONDING  TO  A  GIVEN  TN 


IMPLICIT  REAL*8(A-H,0-Z) 

COMMON/BMATRL/GI , Z , UN , SR , DEBUG 

LOGICAL  FLAG 1 ,FLAG2, DEBUG 

AS-( ( 3 . * ( 1 .-UN)*TNP*SR)/(8.*GI) )*•( 1 ./3.) 

H0»(4.*GI*AS)/t2.-UN) 

RETURN 

END 


SUBROUTINE  ELAD(X,P,DI) 


MATERIAL  DEPENDENT  SUBROUTINE  COMPUTES  THE  ELASTIC 
MODULUS  CORRESPONDING  TO  A  GIVEN  TN 


IMPLICIT  REAL*8(A-H,0-Z> 

COMMON/BMATRL/GI , Z , UN , SR , DEBUG 

LOGICAL  FLAG  1 ,FLAG2, DEBUG 

AS-( (3.*( 1 . -UN)*P*SR)/(8.*GI ))**(!. /3.) 

DI ■ ( 3 .*Z*P*( 2. -UN) * ( 1 .  -( ( 1 .-E/(Z*P) )*( 1 ,-X/(Z*P) ))**(! 
1/3. ) ) )/(8.*GI*AS) 

RETURN 

END 


SUBROUTINE  NORMAD ( EN , TN 1 , TNN , DTN ) 


n 

Lmjiiuiiuuuauiut 


V  % 
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C 

c  COMPUTES  THE  TOTAL  NORMAL  FORCE  AND  THE  PORCE  INCREMENT 

C  BASED  ON  THE  TOTAL  NORMAL  DISPLACEMENT 

C 

IMPLICIT  REAL*8(A-H,0"Z) 

COMOM/BMATRL/GI , Z,UN,SR, DEBOO 
LOGICAL  PLAG 1 , FLAG2 , DEBOG 

TNN-8 . *GI *DSQRT(SR*EM*BH*EN/8 . ) /( 3 . *  < 1 . -UN ) ) 

DTN-TNN-TN1 

RETORN 

END 

C 

C 

C*—*— *********************************** ************************* 

c 

C  SUBROUTINE  CONPLD(M, IAD1 ,TNN,OTN,DEX, DEI, DTE, DTT, IDEE) 

C  THIS  ROUTINE  COMPUTES  THE  INCREMENT  OP  PORCE  POR  A  GIVEN 
C  INCREMENT  OP  DISPLACEMENT  WITHOUT  UPDATING  THE  HISTORY 
C  * 

£•••****••**••****•*••*•*•••**••••**••****••**••••—*•••••••••••••* 

c 

SUBROUTINE  CONPLD ( M , I AO 1 ,TNN , DTN , DEX, DEY , DTE, DTT , IDEE) 

IMPLICIT  REAL*8 ( A-H ,0-Z ) 

COMMON/BMATRL/GI , Z , UN , SR , DEBUG 
COMMON/BHI ST/AHI ST { 1,3004) 

DIMENSION  XXKI  (100),XKRS  (100), HHI  (100), HHS  (100), WX  (100), 
$WY( 1 00) ,TTZS( 1 00) ,TTYS  (100),TTXI  (100),TTYI  ( 100) ,AAXS( 100) , 
SAAYS  (100) ,AAZI  (100),AAYI  (100), IC  (100) 

LOGICAL  PLAG 1 ,PLAG2, DEBUG 
C  CHARACTER* 5  CHECX 

C**************************************************************** 

C - INITIALIZATION - 

C  IP  (CHECK  .NE.  'FIRST* )  THEN 

C  CALL  EMOD(TNN,H0) 

C  CHECK*' FIRST' 

C  END IF 

I DEX-0 
IN-IAD1 
DEBUG-. FALSE. 

NG-IDINT( AHIST(M, IAD1 ) ) 

DTE-0. 

DTY-0. 

NGROUP- 1 00 
C 

C  PREVIOUS  NORMAL  FORCE,  TNI 
C 

TNI -TNN-DTN 
C 

C  MAKE  SURE  THERE  IS  A  CONTACT;  FOR  THIS,  THE  NORMAL  FORCE  HAS  TO 
C  BE  POSITIVE” I DEX- 2  MEANS  THAT  THERE  IS  NO  CONTACT 
C 

IF  (TNN.GT.0)  GO  TO  9827 

I DEX- 2 

RETURN 

9827  IF(NG.NE.O)  GOTO  9820 
C 

C  CALL  THE  ELASTIC  MODULUS  OF  THE  CURRENT  NORMAL  FORCE 
C 


CALL  EMOD(TNN,H0) 


no  non  non  non  non  noon  nnnnjnnn 


tx-ho*dbx 

TT-H0*D«T 

TM-DSQRT ( TX*TX+TT*TT ) 
DTM-TM 

IP(DTN.SQ.O.)  GOTO  9840 
IF(TM.GB.Z*TNN)  THEM 

IDEX-1 

TX-TX*Z*TNN/TM 

TY-TY*Z*TNN/TM 

TM-Z*TNW 

ENOIP 

DTX-TX 

DTT-TT 


INITIALIZE  THE  NEW  CONTACT  HISTORY 


U0  NG-1 


INITIAL  GROUP  CHARACTERISTICS 
- MODULI  AND  SIZES - 


HHS ( 1 ) "HO 
HHI(1 )-0. 

XKKS ( 1 )-Z*TNN 
XKKI ( t >-0. 


-  INITIAL  SURFACE  POSITIONS  - 

COORDINATES  OP  CENTRES  OP  MAX  BOUNDARY  YIELD  CIRCLE 

AAXS(1 )-0. 

AAYS( 1 )"0. 

COORD I ANATES  OP  CENTRES  OF  MIN  BOUNDARY  YIELD  CIRCLE 

AAXId  )-TX 
AAYKO-TY 


INITIALS  REFERENCE  POINTS 


TTXS( 1  )«0. 

TTYS(1 )-0. 

TTXI (1 ) -TX 
TTYId  )-TY 

- KIND  OP  INITIAL  GROUP 

IF(IDEX  .EQ.  0)  THEN 

—  NO  SLIDE  - 

ICC  1 )-0 

wx(  1  )-1  . 

WY  (  1  )-1  . 

VX-1  . 

VY-1  . 

GO  TO  9910 
ELSE 

-  SLIDE  OCCURS  . 
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IC<1)-1 
WX(  t )  -TX/TM 
WT(  1  )-TY/TM 

rrxs( i )-tx 
TTTS( 1 )-TY 
IDEX-1 
GO  TO  9910 
EMOXf 
C 

C -  BEGIN  CALCULATION  WITH  THE  OLD  CONTACT 

C 

9820  DKM-DSQRT(DEX*DEX+DEY*DEY ) 

C 

C—  CHECK  If  THE  MEMORY  IS  ENOUGH - 

C 

I P ( NG . LT . NGROUP )  GO  TO  9790 

IOEZ-4 

RETURN 

C 

C  -  GET  THE  HISTORY  OP  THE  CONTACT  PROM  THE  MEMORY 
C 

9790  CONTINUE 

DDTN-AHIST(M,IAD1+1 ) 

DO  9810  1-1 ,NG 

XAD2-IAD1 +1+1 

XKKI ( I ) -AHI ST (M, I AD2 ) 

KKKS(I ) -AHIST(M, IAD2+NGROUP) 

HHI ( I )-AHIST(M, IAD2+2+NGROUP) 

HHS ( I ) - AHI ST ( M , I AD2+ 3+NGROUP ) 

Will  )-AHIST(M,IA02+4*NGR0UP) 

Will )  -  AHI ST(M,  I AD2+5+NGROUP) 

TTXS ( I ) -AHI ST ( M, I AD2+6+NGROUP ) 

TTYS ( I ) - AHI ST ( M , I AD2 ♦ 7 +NGROUP ) 

TTXI (I)-AHIST(M, I AD2+8+NGROUP) 

TTYI ( I ) -AHI ST (M, I AD2 +9+NGROUP ) 

AAXS ( I ) - AHI ST ( M , I A02 + 1 0 +NGROUP ) 

AAYS ( I ) -AHI ST ( M , I A02+ 1 1 +NGROUP ) 

AAXI ( I ) -AHI ST(M, IAD2+1 2+NGROUP) 

AAYI ( I ) - AHI ST ( M , I AD2+ 1 3* NGROUP ) 

9810  IC(I)-IDINT(AHIST(M,IAD2+14*NGROUP>> 

TX-TTXI  ( NG ) 

TY-TTYI (NG) 


C 

C  - END  OP  INITIAL  REFERENCE  POINTS 

C 

IF(DEM.EQ.0. )  GOTO  9900 
C 

C - UPDATING  PARAMETER  BEFORE  INCREMENT 

C 


HSL-HHS(NG) 

HIL-HHI (NG) 

IP(HIL.LE.0. )  HIL-HSL 
XKL-XKKS(NG) 

I F ( TN 1 . GT . 0 , )  CALL  ELAD ( XKL , TN 1 , EML ) 
HML-0.D0 

I P ( EML . NE . 0 . )  HML-XKL/EML 
NGL-NG 
C 

C -  CALCULATE  THE  CORRECT  FORCE  INCREMENT 

C 


nnn  nnn  «o  nnn  nnnnn 


DTX-HSL*DBX 

DTY-HSL*DRY 

DTM-HSL*DEM 


OMIT  NORMAL  VECTOR 


IP(IC(NOL) . EQ.O)  THEM 

VX-DEX/DEM 

VY-DEY/DSM 

ELSE 

VX-WX(NGL) 

VY-WY(NGL) 

VXL-VX 

VTL-VT 

ENOIP 

DBN-DTX*VX>DTT«VT 
DBT-DABS ( DTX*VY-DTY*VX ) 

PINO  THE  CORRESPONDING  FORCE 

I P ( DBN . LE . Z*DTN )  GOTO  9900 
IFdC(NGL).EQ.O)  THEN 
DTM-(DEM*Z*DTN*(1 ./HIL-1 ./HSL) )*HIL 
ELSE  IP(HML.EQ.O.)  THEN 

DTM-DEM* ( DEM* Z  *DTN* ( 1 . /H I L- 1 . /HSL ) ) *H I L/ ( DEX*VX*DEY*VY ) 
ELSE 

DTM-DEM* ( DEM* Z*DTN* ( 1 ./HIL-1 . /HSL ) )/( ( DEX*VX+DEY*VY ) /HI L 
$  +DABS(DEX*VY-DEY*VX)/HML) 

ENOIP 

OTX-OEX*OABS ( DTK) /DEM 
DTY -DEY*OABS ( DTM ) /DEM 
DBN-DTX*VX+DTY*VY 
DBT-DABS ( DTX*VY-DTY*VX ) 

10  TN-TNN 
TX-TX+DTX 
TY-TY+DTY 

TM-DSQRT (TX*TX+TY*TY ) 


*  -■*  j- 

■V 

.V  A 

■V 


- UPDATING  FORCE  INCREMENT 

TNL-TN 
TXL-TX 
TYL-TY 

- CHECK  IP  THERE  IS  FAILURE 

IF(TM.LT.Z*TN)  GO  TO  9910 
IDEX-1 

TX-TX*Z*TN/TM 
TY-TY*Z*TN/TM 
TM-Z*TN 

DTX-TX-TTXI (NG) 

DTY-TY-TTYI ( NG ) 

GOTO  9840 
C 

9910  CONTINUE 
RETURN 
END 


wwpaygwwwtwiwmTOn  x."  x_t^  it"  ^ v  W '-  Tv  w  v  ▼' 


v'-  v\  t;v;^.^’V'  W  wJTWk¥ »  }T*  V 
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c 

c 

C  THIS  SUBROUTINE  COMPOTES  THE  INCREMENT  OF  DISPLACEMENT  FOR 
C  A  GIVEN  INCREMENT  OF  FORCE  WITHOUT  STORING  OR  UPDATING  THE 
C  STRESS  HI STORE. 

C 

C***************************************************************** 

c 

SUBROUTINE  CONPLF(M, IAD! ,TNN,DTN, DTE, DTY,DBX,DBY, IDEE) 
IMPLICIT  RSAL*8(A-H,0-X) 

DIMENSION  XKXI  (100)  ,XKKS  (100),  HHI  (  100)  ,HHS(  100)  ,WX(  100) , 
$WT  (tOO) ,TTZS( 100) ,TTTS( 100) ,TTZZ (100) ,TTTI  { 100) ,AAXS( 100) , 
SAATS  (100) , AAAI  (100),AATI  (100), IC  (100) 

COM40N  /BMATRL/  GI , Z ,UN, SR, DEBUG 
COHMON/BHI ST/  AHIST< 1 , 3004 ) 

LOGICAL  FLAG 1 .FLAG 2, FLAG 3, DEBUG 


C  CHARACTER* 5  CHECK 

c - INITIALISATION 

C  IF  (CHECK  .NE.  ’FIRST*)  THEN 

C  CALL  EMOD(TNN,H0 ) 

C  CHECK- 'riRST' 

C  END IF 


NG-IDINT( AHIST(M, I ADI ) ) 

NGROUP- 1 00 
IDEX-0 
C 

C  PREVIOUS  NORMAL  FORCE, TNI 
C 

TN 1  -TNN-DTH 

C -  MAKE  SURE  THAT  THE  NORMAL  FORCE  IS  POSITIVE.  OTHERWISE, 

C -  IDEX-2  AMD  THERE  IS  NO  CONTACT  BETWEEN  THESE  PARTICLES 

I F ( TNN . GT . 0 . )  GOTO  9827 

I DEE- 2 

RETURN 

9827  IF(NG.NE.O)  GOTO  9820 

C  CALCULATE  THE  ELASTIC  MODULUS  OF  THE  CURRENT  NORMAL  FORCE 
CALL  EMOO(TNN,HO) 

TX-DTX 

TY-DTY 

TM-DSQRT ( TX*TX*TY*TY ) 

DTM-TM 

IFtDTM.EQ.O. )  GOTO  9840 


C -  IF  THE  FORCE  POINT  IS  OUTSIDE  THE  FAILURE  SURFACE  - 

C - USE  THE  FORCE  POINT  ON  THE  FAILURE  SURFACE . . - 


IF(TM.GE.Z*TNN)  THEN 
IDEX-1 

TX«TX*TNN/TM 

TY«TY*TNN/TM 

TM-Z*TNN 

END  IF 

DTX-TX 

DTY-TY 

DEX-TX/H0 

DEY-TY/H0 


C -  INITIALIZE  THE  NEW  CONTACT  HISTORY 

9840  NG-1 

C -  INITIAL  GROUP  MODULI  AND  SIZES  — 


UU  o  O  u  U  U  UOWUUO  co  uou 
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HHKO-O. 

XMSd  )«Z*TNN 
XRRKD-O. 

-  INITIAL  YIELD  SURFACE  CENTERS  - 

COORDINATES  OP  CENTRES  OP  MAX  BOUNDARY  YIELD  CIRCLE 
AAXS(1)-0. 

AAYS(1 )-0. 

COORD I ANATSS  OP  CENTRES  OP  MIN  BOUNDARY  YIELD  CIRCLE 
AAXId  )-TX 
AAYId)«TY 

- INITIAL  REPERENCS  POINTS  - 

TTXS(t)-0. 

TTY5(1 )-0. 

TTXI (1 ) -TX 
TTY  I  (  D-TY 

- kind  op  initial  group  - 

IP(IDEX  .EQ.  0)  THEN 

-  NO  SLIDE  - 

IC(1 )-0 

wxd)-i. 

WY(  1  )-1  . 

VX-1  . 

VY-1  . 

GO  TO  9910 
ELSE 

-  SLIDE  OCCURS  - 

IC(1 )-1 
WX(  1  )  -TX/TM 
WY  (  1  )  -TY/TM 
TTXS ( 1 ) “TX 
TTYSd  )-TT 
IDEX-1 
GO  TO  9910 
ENDIP 

-  CALCULATION  BEGIN  WITH  OLD  CONTACT  - 

9020  DEM-DSQRT(DEX*DEX+DEY*DEY) 

020  CONTINUE 

-  STOP  WHEN  THE  MEMORY  IS  NOT  ENOUGH  - 

IP(NG.LT.NGROUP)  GO  TO  9709 
WRITE ( 3,0101 )  NG 

101  FORMAT ( 1  THE  MEMORY  OP  THE  YIELD  SURFACE  IS  NOT  ENOUGH' ,14) 
STOP 

GET  THE  CONTACT  HISTORY  FROM  MEMORY 

9709  CONTINUE 

DDTN ■ AH I  ST ( M ,  I  AD  1  + 1  ) 

DO  9010  1-1 ,NG 

IAD2-IAD1 ♦ 1 

XKKI ( I ) -AHI ST(M, I AD2 ) 

XKKS ( I ) - AHI ST ( M , I AD2+NGROUP ) 

HHI (I )«AHIST(M, I AD2+2*NGROUP ) 

HHS ( I ) -AHIST(M, IAD2+3*NGROUP) 

WX(I )«AHIST(M,IAD2*4*NGROUP) 

WY  (I ) -AHI ST(M, I AD2  +  5*NGROUP ) 

TTXS( I ) -AHIST(M, I AD2+6*NGROUP) 

TTYS(I ) -AHI ST(M, I AD2+7*NGROUP) 
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TTXI (I ) -AHISTtH, IA02+8*NGR0UP) 

TTTI ( I ) -AHI ST(M, I AD2*9*NGROUP) 

AAXS ( I ) -AHISTtM, IAD2+1 0*NGROUP) 

AATS(I )«AHIST(M, IAD2+1 1 *NGROUP) 

AAXX d )-AHIST(M, IAD2+1 2*NGR0UP) 

AAYI (I )“AHIST(M, IAD2+1 3*NGR0UP) 

IC(I )-IDINT(AHIST(M, IAD2-M 4*NGR0UP) ) 

9810  COHTI HOT 

TX-TTXI(NG) 

TT-TTTX(HG) 

C 

C - UPDATI MG  PARAMETER  BEFORE  INCREMENT - 

C 

HSL-HHS(NG) 

HXL-HHI(NG) 

I F ( HI L . LE . 0 . )  HIL-HSL 
XKL-XKKS(NG) 

I F ( TN 1 . GT . 0 . )  CALL  ELAD ( XXL , TN 1 , EML ) 

HML-0.D0 

I F ( EML . NE . 0 . )  HML-XXL/EML 
NGL-NG 

C -  CALCULATE  THE  CORRECT  FORCE  INCREMENT  - 

DTM-DSQRT ( DTX*DTX*DTY*DTY ) 

I F ( DTM . EQ . 0 . )  GOTO  9900 

C - UNIT  NORMAL  VECTOR - 

C*********** ********************** ********************** 
IF(IC(NGL) .EQ.O)  THEN 
VX-DTX/DTM 
VT-DTT/DTM 
ELSE 

VX-WX(NGL) 

VY-WT(NGL) 

VXL-VX 
VTL-VT 
END  IF 

DBN-DTX*VX*DTY*VT 
DBT-DABS ( DTX*VY-DTT*VX ) 

C -  TO  FIND  THE  CORRESPONDING  DISPLACEMENT  - 

IFfDBN.LE. Z*DTN)  THEN 

DEX-DTX/HSL 

OEY-DTY/HSL 

ELSE  IF(IC(NGL). EQ.O. OR. HML. EQ.O.)  THEN 
DEX- Z *DTN*VX/HSL* ( DBN- Z*DTN ) *VX/HI L 
DEY- Z  *DTN*VY/HSL+ ( DBN- Z  *DTN ) • VY/H I L 
ELSE 

DEX»Z*DTN*VX/HSL«-(DBN-Z*DTN)*VX/HIL+DBT*VX/HML 
C  DEX-Z*DTN*VX/HSL*(DBN-Z*DTN)*VX/HIL+(DTX-DBN*VX)/HML 
DE7- Z  *DTN*VY/HSL+ ( DBN- Z  *DTN ) *V7/HI L* ( DTY-DBN*VY ) /HML 
END  IF 

DEM-DSQRT( DEX*DEX*DEY*DEY ) 

GOTO  9900 

C -  CALCULATION  OF  THE  CURRENT  YIELDING  SURFACE  BEGIN  --- 

9900  TN-TN1+DTN 
TX-TX+DTX 
TY-TY*DTY 

TM«DSQRT( TX*TX+TY*TY ) 

C - UPDATING  FORCE  INCREMENT - 

TNL-TN 

TXL-TX 


CHECK  IF  THERE  IS  FAILURE' 


IF(TM.LT.Z*TN)  GO  TO  9910 

IDEX-1 

TX«TX*Z*TN/TM 

TT-Tr*X*TN/TM 

TM-X*T* 

DTX-TX-TTXI (NO) 

DTT-TT-TTTI (HO) 

DBM»DTX*VX+DTT*VY 

DBX-X*DTH*VX/HSL* ( DBN- Z*DTN ) *VX/HSL 
DET-X*DTN*VT/HSL* ( DBN- Z*DTN ) •VT/HSL 
IF(IC(NGL)  .EQ.  0  .OR.  HML  .EQ.  0.)  THEN 
GO  TO  9840 
ELSE 

DEX-DEX* ( DTX-DBN*VX ) /HML 
DEY-DET-*- ( DTT-DBN*VT ) /HML 
ENDIF 

GO  TO  9840 
9910  CONTINUE 
RETURN 
END 


APPENDIX  D 

Media  Configurations  Used  to  Simulate 
the  Aggregate 
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Figure  Dl.  Configuration  of  the  Media  Used  in  the  Simulations 
and  Orientation  of  Each  Element  . 
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Figure  D2.  Configuration  of  the  Media  Used  in  the  Simulations 
and  Orientation  of  Each  Element  . 
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Figure  D3.  Configuration  of  Medium  5  Used  in  the  Simulations 
and  Orientation  of  Each  Element  . 


P 

£ 


ft 

*— ■ 

8 

a 


immmm 


Medium  2 


APPENDIX  E 

Compression  with  Constant  Mean  Stress. 
Element  Stress-Strain  Behavior 


220 


MEW 


— '  0  32 


0.00  - 1 - 1 - 1 _ I _ 1 _ | _ I _ L-J 

0.00  0.04  0.08  0.12  0.18  0  20 

1 60*-  y 1 2  (x  10‘3) 


128 

(M 

E 

P  0.96 


b  0  32 


0  00*— ^ - 1 - 1 - 1 1— J  1111 

000  0.08  0.18  0.24  032  0.40 

1 60r  £22  (x  10“3) 


E 

O  0.96 
« 

O) 

■X  0.64 


b  0  32 


0  00  0.08  0.18  0  24  0  32  0  40 


fit  (x  icr3) 


Figure  Eli.  Medium  2:  Compression  with  Constant  Mean  Stress 
(oq*1.0  Kg*/ cur).  Stress-Strain  Behavior  of 
Element  Ell  Oriented  at  8*20°. 
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Figure  E14.  Medium  2:  Compression  with  Constant  Mean  Stress 
(ao“i.O  Kg*/cnr).  Stress-Strain  Behavior  of 
Element  E14  Oriented  at  8*60°. 
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APPENDIX  F 


Biaxial  Compression  with  Variable  Mean  Stress 
Element  Stress-Strain  Behavior 
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Figure  F5.  Medium  1:  Compression  with  Variable  Mean  Stress 
(oq-I.O  KgVcm*).  Stress-Strain  Behavior  of 
Element  5  Oriented  at  S“80°. 
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Figure  F9.  Medium  1:  Compression  with  Variable  Mean  Stress 
(oq*1.0  Kg*/cm^).  Stress-Strain  Behavior  of 
Element  9  Oriented  at  Sa20°. 
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Figure  F13.  Medium  1:  Compression  with  Variable  Mean  Stress 
(ag»1.0  l^*/cm‘).  Stress-Strain  Behavior  of 
Element  13  Oriented  at  6*60°  . 
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Medium  1:  Compression  with  Variable  Mean  St 
(oq“1.0  Kg*/cm2).  Stress-Strain  Behavior  of 
Element  14  Oriented  at  0*50°  . 
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APPENDIX  G 


Medium  2:  Compression  with  Mean  Stress  Constant 
Variation  of  Normalized  Stress  and 
Strain  During  Loading 
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Figure  G2.  Medium  2:  Compression  with  Constant  Mean  Stress 
Oq«1.0  Kg*/cm2,  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  20°. 
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Figure  G5. 


Medium  2:  Compression  with  Constant  Mean  Stress 
Oq-1.0  Kg*/cm2.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  50s. 


sM? 


fi-fa(x  10"3) 


0.00  0.04  0.08  0  12  0.16  0  20 

i  oop  fi-ea(x  tO'3} 


0.00  0.04  0.08  0.12  0.16  0  20 

fi-£a(x  10'3) 


Figure  G6, 


Medium  2:  Compression  with  Constant  Mean  Stress 
Oq«1.0  KgVcrn^.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  60°. 
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Figure  G9. 


Medium  2:  Compression  with  Constant  Mean  Stress 
Oq»1.0  KgVcm^.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  10°. 
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Figure  Gil.  Medium  2:  Compression  with  Constant  Mean  Stress 
Oq«1.0  f^j*/cm^.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  30°. 
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Figure  G12.  Medium  2:  Compression  with  Constant  Mean  Stress 
ao«1.0  Kg*/cm^.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  40®. 
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Figure  G13.  Medium  2:  Compression  with  Constant  Mean  Stress 
ao*1.0  Kgtycm^.  Normalized  Element  Stress  Versus 
Applied  Principal  Strain  Difference  for  All 
Elements  Oriented  50°. 
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